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nozzle cross-sectional area, m 2 

constant coefficients in equation (36) for C p ^ 

initial approximation of mols of CO 2 per total mol equiva- 
lence of C in mixture 

constant coefficients in equation (35) for K p j 
molar heat capacity of species k } J/mol-K 
specific heat capacity of carrier gas, J/kg-K 
specific heat capacity of species fc, J/kg-K 
constant coefficients in equation (63) for 
constant coefficients in equation (66) for A ^ 
function defined by equation (111) 

factor to account for carrier gas defined by equation (74) 

function defined by equation (46) 

specific enthalpy of total mixture, J/kg 

specific enthalpy of gaseous mixture, J/kg 

specific enthalpy of liquid water, J/kg 

stagnation specific enthalpy, J/kg 

specific enthalpy of fuel, O 2 , and N 2 in feed, J/kg 

total specific enthalpy of feed to combustion chamber, J/kg 

specific enthalpy of evaporation at 373.15 K, J/kg 

nucleation rate, droplets formed/m 3 -s 

mean nucleation rate defined by equation (108), droplets 
formed/m 3 -s 

equilibrium constant for reaction j 
Knudsen number defined in terms of r 
Knudsen number defined in terms of r* 

Boltzmann constant, J/molecule-K 
latent heat of evaporation at T 5 (pi), J/kg 
mean free path, m 
mass of a water droplet, kg 

average molecular mass of gaseous mixture, kg/molecule 
average molecular mass of carrier gas, kg/molecule 
molecular mass of water, kg/molecule 
mass flow rate, kg total mixture/s 

iteration index and atoms of carbon per molecule of fuel 
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t g 

Tl 

Tr 

T c 
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number of droplets of kind i per mass of total mixture, 
droplets/kg 

Prandtl number defined in terms of c Pj i (eq. (69)) 
pressure, N/m 2 

partial pressure of water vapor, N/m 2 

vapor pressure of water over a flat surface, N/m 2 

standard state pressure, 0.101325 MN/m 2 

Kantrowitz correction factor, equation (62) 

condensation coefficient 

evaporation coefficient 

universal gas constant, 8.314 J/mol-K 

specific gas constant for mixture, R/W , J/kg-K 

universal gas constant, 1.987 cal/mol-K 

specific gas constant for carrier gas, R/W Ci J/kg-K 

specific gas constant for water, R/W\, J/kg-K 

droplet radius, m 

nozzle wall radius, m 

critical droplet radius, m 

specific entropy of total mixture, J/kg-K 

specific entropy of water vapor at standard state pressure, 
J/kg-K 

specific entropy of liquid water, J/kg-K 

specific entropy of evaporation at standard state pressure, 
J/kg-K 

temperature, K 

temperature of gaseous mixture, K 

temperature of liquid water, K 

reduced temperature, T jT c 

critical temperature of water, 647.3 K 

saturation temperature based on partial pressure of water 
vapor, K 

time, s 

time step, s 

velocity, m/s 

average molecular weight of gaseous mixture, kg/mol 
(also W — 1/Y) 

average molecular weight of carrier gas, kg/mol 
molecular weight of species fc, kg/mol 
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Subscripts: 

c 


mass of liquid water per mass of total mixture 
mass of carrier gas per mass of total mixture 
distance along nozzle measured from throat, m 
computational step along nozzle, m 
mol number of mixture, mols/kg 

equivalent total mols of elements C, H, O, and N per mass 
of total mixture, mols/kg 

mol number of species fc, mols/kg 

mol number of species k in gaseous mixture at point in 
nozzle where J = J m j n , mols/kg 

mols of fuel, O 2 , and N 2 per mass of total mixture in feed, 
mols/kg 

mol fraction of species k 
parameter defined by equation (112) 
ratio of droplet radius to critical radius, r/r* 
constant in equation (78) 

thermal accommodation coefficient for carrier gas 
interaction with water droplets 

constant in equation (75); also, Langmuir parameter 

ratio of specific heats 

parameter defined by equation (82) 

nozzle boundary layer displacement thickness, m 

allowable iteration error 

elemental ratio of hydrogen to carbon in mixture 

elemental ratio of nitrogen to oxygen in mixture 

parameter defined by equation (75) 

parameter defined by equation (76) 

thermal conductivity of gaseous mixture, J/s-m-K 

viscosity of gaseous mixture, N-s/m 2 

parameter defined by equation (78) 

parameter defined by equation (77) 

mass density, kg/m 3 

surface tension, N/m 

equivalence ratio 

parameter defined by equation (80) 
parameter defined by equation (73) 

carrier gas 
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gaseous mixture 

properties at interface between free molecular and contin- 
uum regimes (see appendix A) 

location index along nozzle axis 

chemical species index (where k = 1 denotes water vapor) 

liquid water 

minimum 

summation index 

stagnation conditions 

saturation condition except in equation (12) where it is 
stoichiometric condition 

water vapor 

evaluated at 373.15 K 


Superscripts: 

i droplet type index 

o conditions at point in nozzle where nucleation just 

begins, except p° denotes standard state pressure 
of 0.101325 MN/m 2 

f approximation 

ff alternate approximation 

* sonic condition at nozzle throat 

A caret (A) over a symbol indicates the feed condition to combustor. 
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Introduction 


Experimental studies of propulsion systems and 
aerothermal-structural systems concepts for hyper- 
sonic flight require wind tunnels that generate high 
enthalpy flow. One approach to achieving such 
flows utilizes combustion at a relatively high pres- 
sure followed by expansion of the resulting combus- 
tion products to form the test stream. A number of 
combustion-driven wind tunnels burning various fu- 
els and of various sizes have been built and operated. 

Fuels for these facilities include hydrogen, 
methane (as the major component of natural gas), 
and isobutane each burning with air or oxygen- 
enriched air. In all cases, the products of combustion 
contain a substantial amount of water vapor. Con- 
ditions for which condensation of water could occur 
depend on the equivalence ratio, temperature, and 
pressure in the tunnel combustion chamber, and the 
extent of the expansion of the combustion products 
in the nozzle. If liquid water does form, the process 
experiences entropy production and the flow proper- 
ties are altered relative to an isentropic expansion 
There is therefore a need to be able to predict and 
analyze the effects of nucleation and water droplet 
growth in combustion-driven wind tunnels. 

A great deal of computational and experimental 
work has been carried out over the years on the nu- 
cleation and droplet growth of liquid water from pure 
steam. References 1 and 2 include reviews of previ- 
ous studies and present computations and compar- 
isons with experimental data for pure steam. One 
of the most recent studies has been presented by 
Young in references 3 and 4 and is the starting point 
and basis for the present work. The method pre- 
sented by Young for steam is modified and extended 
to the rapid expansion of combustion products. This 
requires the addition of a scheme to compute the 
properties of the reacting gas mixture in the com- 
bustion chamber and the subsequent expansion along 
the nozzle prior to the formation of liquid water and 
the modification of the nucleation and droplet growth 
equations to account for the noncondensable compo- 
nents of the combustion products. 

The equations are developed for a quasi-one- 
dimensional flow of combustion products; nucleation 
and droplet growth are taken into account in the de- 
velopment. A numerical scheme that employs these 
equations is presented. Sample results for selected 
conditions in the Langley 8-Foot High- Temperature 
Tunnel are also given and show the general effects of 
water condensation. 


Development of Equations 

The development of an appropriate set of equa- 
tions for water condensation from a rapidly expand- 
ing flow of combustion products follows closely the 
approach presented by Young (refs. 3 and 4) for con- 
densation of pure steam. In the present analysis, the 
tunnel stagnation conditions and the subsequent ex- 
pansion along the nozzle are described by a reacting 
mixture of ideal gases in chemical equilibrium. The 
condition of chemical equilibrium is applied along the 
expansion until nucleation just begins, at which point 
the temperature is low enough so that no further re- 
action needs to be considered. Also the nucleation 
and droplet growth processes take place in a gaseous 
mixture of water vapor and a noncondensable car- 
rier gas composed of the other combustion products 
rather than pure steam. It should be noted that no 
liquid water is formed while the condition of chemical 
equilibrium is imposed and that no chemical reaction 
takes place during nucleation and droplet growth. 
Figure 1 is a schematic of a combustion-heated wind 
tunnel to which this analysis is applied. 

The description of the nozzle flow process from 
the combustion chamber through the throat and un- 
til nucleation just begins requires a set of flow equa- 
tions, stoichiometric relationships, a set of equations 
for computing the equilibrium chemical composition, 
thermodynamic relations and properties for a react- 
ing gas mixture, and a nucleation rate equation. This 
process is indicated as part 1 and part 2 in figure 1. 
Computation of the flow process for continued expan- 
sion during nucleation and droplet growth requires 
additional equations for transport properties, droplet 
growth rate, and entropy production. This process 
is indicated as part 3 in figure 1. The following 
paragraphs present the specific equations required for 
these processes. 

Flow Equations 

The continuity equation for quasi-one-dimensional 
nozzle flow is 


. PgAU 

m = 

1 — w 

where it has been assumed that liquid water and the 
gaseous mixture have the same velocity. Equation (1) 
applies along the entire expansion prior to liquid 
water formation with w = 0 as well as with liquid 
water present. 

The momentum equation is 


(l-w)d P = _ udu 

PG 


( 2 ) 



where again the liquid water and gaseous mixture 
have the same velocity. 

The energy equation for adiabatic steady flow is 

h 0 = (l-w)h G + J2 wihi L + ^- (3) 

i 

where the first term on the right-hand side of equa- 
tion (3) can be written as 


(1 - w)h G - (1 - w - w c )h G l + w c h G)C (4) 

Now in the region where water vapor condenses and 
no further chemical reaction occurs in the gas phase, 
the differential form of the energy equation (eq. (3)) 
becomes 


(1 — w — w c ) dhQ^i 4- w c dliQ^ -f U dU 

= h G ,i dw — w l dh l L - h l L dw l (5) 

i i 


Stoichiometric Relationships 

The computation of the chemical composition of 
a reacting gas mixture composed of the elements C, 
H, O, and N requires a specification of the elemental 
constants Yq, Yp, Yq , and Y n which represent the 
equivalent total number of mols of each element 
per mass of mixture. These four constants can be 
determined by specifying the elemental nitrogen-to- 
oxygen ratio the elemental hydrogen-to-carbon 

ratio ?7 h,C> the equivalence ratio <£, where 

( 6 ) 

and the four elemental constants are related by the 
identity 


12Y C + Y H + 16Yo + 1 4F N = 1000 (7) 

It follows then that 


Y 0 = 


500 


8 + 7? ?N,0 + 4 >\{ 12 + f?H,c)/( 4 + ^H.c)] 
24>Yo 


4 + ^?H,C 
= ^N.O^O 

Y H = f?H,C*C 


(8) 

(9) 

( 10 ) 

(11) 


These elemental constants will be used in the set of 
equations appearing in the next section. 

The chemical composition of the feed stream to 
the combustion chamber is also required. Although it 
is possible to present a scheme for a more general fuel 
type, attention is limited to aliphatic hydrocarbons 
burning with a mixture of oxygen and nitrogen. For 
this class of fuels, the stoichiometric reaction is 

C„H m + (n + j) 0 2 = nC0 2 + j H 2 0 


and the stoichiometric ratio of fuel to oxygen in the 
feed is 


(SY | = i 

\Y 0 J a ^ + m,c) 

so that 


Yf 


40 

n (4 + ^7H,c) 


?o 2 


Also of course 


( 12 ) 


(13) 


^N 2 = W,0*b 2 (14) 

and since 

?O a - ^ (15) 

it follows that 


° 2 8 + 7»7n,o + <t> [(12 + ^H.c)/^ + rm,c)} 

(16) 

In the flow region where nucleation and droplet 
growth occur, no further chemical reaction takes 
place so that the mixture composition changes only 
as a result of water condensation. The gaseous 
mixture composition in this region can be expressed 
in the following way in terms of the mass fraction 
of liquid water formed and the mol numbers of each 
species in the gaseous mixture just before any liquid 
water is formed. Let w represent the mass of liquid 
water in the stream per unit mass of total mixture 
(liquid water plus gaseous mixture). Also let Y 
represent the mol number of species k in the gaseous 
mixture and Y£ denote the mol number of species k 
in the gaseous mixture just before any liquid water is 
formed. Now since the number of mols of liquid water 
formed per unit mass of the total mixture is w/W\ 
and the mass of gaseous mixture per mass of total 


2 



mixture is (1 — it;), it follows that the mol number of 
water vapor in the gaseous mixture is 


Yi = 


Yi - wjW\ 
1 — w 


(17) 


H 2 , N 2 , H, O, OH, and NO. The composition is de- 
termined by the simultaneous solution of the follow- 
ing six chemical equilibrium relations corresponding 
to the six listed reactions: 

For reaction I, CO 2 + H 2 = CO + H 2 O, 


in the region where nucleation and droplet growth 
occur. The mol numbers for the other species, that 
is for k > 1, are 


*p,i = 


YxY 3 


_ Y k° 


y * = T 


W 


It also follows, since 


W = 


1 


10 

E Y k 

k = 1 


(18) 


(19) 


Y 2 Y 5 

For reaction II, 2 CO 2 = 2CO + O 2 , 


K 


YjY4 P 

p ’ 2 y 2 y p° 


For reaction III, H 2 + O 2 = 20H, 

K. — 


3 


that the average molecular weight of the gaseous 


T4I5 


(25) 


(26) 


(27) 


mixture is 


For reaction IV, H 2 = 2H, 


£ 
I 1 

0 

11 

(20) 

K - y 8 P 

p ’ 4 y 4 y p° 

(28) 

where 




10 

y° = Y, n° 

k= 1 


For reaction V, O 2 = 20, 


(21) 

K Y i P 

p>5 y 4 y p° 

(29) 

Also, the mol fraction of any species in 
mixture is 

the gaseous 

For reaction VI, O 2 + N 2 = 2NO, 


Vk = Y k w 

(22) 

«s £ 
11 

CD 

(30) 

The average molecular weight of the carrier gas, 
that is for all species in the gaseous mixture except 

and the four elemental balance equations: 


water vapor, is 




T „ ^ Y k W k 

w -EyHi 

(23) 

y H = 2 y 4 4 - 2 y 5 + y 7 + y 9 

(31) 

k= 2 


y 0 = y 1 + 2 y 2 + y 3 + 2 y 4 + y 8 + y 9 + y 10 ( 32 ) 

or noting that Y = 1 jW and using the foregoing 
expressions for i/*., Y\ , Y*., and IT, it follows that 

y N = 2 y 6 + y 10 

(33) 

Hr - T y °* Wk 

c ~h ‘-t 

(24) 

y c = y 2 + y 3 

(34) 

which is of course independent of the amount of water 

The expression for K p j is 



condensed, w. 

Equilibrium Chemical Composition 

The equilibrium chemical composition of combus- 
tion products includes 10 species numbered from 1 
to 10 in the following order — H 2 O, CO 2 , CO, O 2 , 


tfpj = (e B ".J r S' 2 ) ( 35 > 

where the constant coefficients B n j are given in 
appendix B. The constants Yq, Yh> an< ^ ' m 
the four elemental balance equations are determined 


3 



illy ill m mmm m i h mi mu mi 


from equations (8) through (11) and the specification 
of *7N,Cb VU } Ci and 0* 


Thermodynamic Relations and Properties 

The starting point for determining the thermo- 
dynamic properties of a mixture is an equation for 
the molar heat capacity at constant pressure of each 
chemical species expressed as 


6 

Cp,k = £ An,kTZ~ 2 (36) 

n= 1 

where A n is a set of constant coefficients given 
in appendix B for each species over a temperature 
range. In some cases, the specific heat capacity at 
constant pressure is useful and can be written as 


c p,k ~ 


w k 


(37) 


In addition to this expression for the kth species, 
expressions for the specific heat capacity at constant 
pressure and the ratio of specific heats for water 
vapor are needed; that is, 


and 


c p,l — 


E%=i A nA T ( 

W i 


n — 2 


71 = 


c p, 1 


c p,l ^1 


(38) 

(39) 


Similar expressions for the carrier gas are also needed 
for computations of droplet growth and can be writ- 
ten as 


10 

E y° k c p , k 


— k— 2 

p ’ c (1 - y°)W c 

(40) 

_ C P.« 

7 c — -p 

C p,c 

(41) 


The specific enthalpy of the gaseous mixture can 
be obtained from the equation for C p integration, 
and summation over all species to give 


10 


h G = J2 Y k 


k= 1 


Ai,k In Tg + V' G . + A 7 

Z — ^ n — 1 


n=2 


(42) 


where the additional constant coefficients in- 

clude the standard state enthalpy of formation of 
each species at the reference temperature of 298.15 K. 


Numerical values of Aj ^ are given in appendix B. 
The units for Hq are joules per kilogram of gaseous 
mixture. An expression for the specific enthalpy of 
water vapor in the mixture is also needed; that is, 


h G , l - ^ 


^1,1 In T G + £ 

n—2 

wT 


r»Tl — 1 


^ n , 1 1 G 


n i 




(43) 


The units of \ are joules per kilogram of water 
vapor. 

The specific entropy of the gaseous mixture can 
also be obtained from the expression for C p inte- 
gration, and summation over all species to give 


10 

=Y, Y k 

fc= 1 


r -a 


1 ,k 


T G 


+ (A 2>k - R') In T G 


6 A rpTl — 2 


n=3 


~R In P G 

_ /D\ 

- ^ ln {^)~ R E Y k 1“ Y k (44) 


where Ag^ includes the standard state entropy of 
formation of species k at the reference temperature of 
298.15 K. Numerical values are given in appendix B. 
For some applications it is convenient to invert this 
equation to express pq as a function of sq , Tq , and 
the set of mol numbers Y^\ that is, 


PG = exp 


G - s G - R In (g/p°) - R t Yk In 
_ 


(45) 

where 


10 

G=E y k 

k= 1 

6 


l l,* 


T G 


+ { A 2 ,k - R') In Tq 


^ A n k TVT 2 


n= 3 


(46) 


The equation of state for the gaseous mixture is 


P = p g RT g 


(47) 


where R = and W — [ Y1 


10 

E 

ik= 1 


-1 
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Computations in the flow region where nucle- 
ation and droplet growth can occur require ad- 
ditional equations for the enthalpy and entropy 
of liquid water. An equation for the specific 
enthalpy of liquid water is obtained by starting 
with the specific enthalpy of evaporation of wa- 
ter at 1 atm and 373.15 K from reference 5; 
that is, A h° L G 373 15 = 2.2570 x 10 6 J/kg. Then us- 
ing equation (43) at T = 373.15 K and the coef- 
ficients A n> i given in table Bl(a), it follows that 

^<7,1,373.15 “ —13.2933 x 10 6 J/kg. Also, the heat 
capacity of liquid water in the region of interest is es- 
sentially constant and is C£ = 4.2 x 10 3 J/kg-K. Now 
since 

h L - h Lt 373.15 = (4.2 x 10 3 )(r - 373.15) (48) 

and 

& h L,G, 373.15 = h G,l, 373.15 ~ h L, 373.15 (49) 

it follows that the specific enthalpy of liquid water is 


where the units of s l are joules per kilogram of liquid 
water-kelvin. 

An equation for the specific latent heat of evap- 
oration of water at a temperature T is also needed 
and is 

L = h G> i - h L (55) 

It then follows from equation (50) that 

L = h G>1 -4.2 x 10 3 r+ 17.1175 x 10 6 (56) 

where hQ\ is evaluated at T and the units of L are 
joules per kilogram of water. 

The vapor pressure of water is also required. An 
expression for the vapor pressure of water over a flat 
surface as a function of temperature was obtained 
by fitting a curve to the tabulated values given in 
reference 5; that is, 

( 664 1 7 \ 

55.897 --4.4864 In T j (57) 


h L = (4.2 x 10 3 )T - (17.1175 x 10 6 ) (50) 

where the units of h l are joules per kilogram of liquid 
water. 

The specific entropy of liquid water can be ob- 
tained in a similar way. The specific entropy of evap- 
oration of water at 1 atm and 373.15 K is 


A4.G.373.15 = A/t 373.i5 3 1 ~ = 6 0485 x 103 J / k S' K 

The entropy of water vapor at 1 atm and 373.15 K 
is obtained by using the equation 


-*Ai,i/Tg + ^ 2,1 In T g + n 'n-% f Ag t i 

o . n=3 

s c.i = 

(51) 

and the coefficients A n? i given in table Bl(a) so that 
s°q l 373 15 = 10.905 x 10 3 J/kg-K. Now since 

S L ~ 5 L, 373.15 = ( 4 - 2 X 1° 3 ) l n ^ ( 52 ) 

and 

G, 373.15 = S G, 1,373.15 ” 5 L, 373.15 ( 53 ) 

it follows that the specific entropy of liquid water is 

s L = 4.2 x 10 3 In T - 20.0149 x 10 3 (54) 


where p 0 0 has units of newtons per meter 2 . 

The surface tension of water is represented by the 
equation given by Young (ref. 4); that is, 

cr = (82.27 + 75.612T# - 256.889T^ 

+ 95.9287^) x 10" 3 (58) 

where 7/2 = ^- with T c denoting the critical tem- 
perature of water. The units of o are newtons per 
meter. 

Finally, the critical droplet radius is obtained 
from the well-known equation 


PlRiT g In ( pi/poo ) 

where p\ = y\p and r* has units of meters. Droplets 
larger than r* tend to grow while droplets less than 
r* tend to evaporate. 

Nucleation Rate 

The nucleation rate equation given by Young can 
be modified to take into account a carrier gas by 
replacing the gas density in Young’s equation which 
is for pure steam by the density of water vapor. This 
follows from the argument that nucleation depends 
on collisions between water molecules in the gaseous 
phase and clusters of water molecules, both of which 
are proportional to the density of water vapor. Now 
since the density of water vapor in a gaseous mixture 
can be expressed as 
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PG, 1 


vim 

m 


PG 


(60) 


the nucleation rate equation for a gaseous mixture 
that includes a carrier gas becomes 


The thermal conductivity of the gaseous mixture 
can also be obtained by starting with a polynomial 
expression in terms of temperature for each pure 
species 


J = 


<? c 

/yi^n 2 1 

{ 2cr \ 1//2 Pq 

1 + Q 

(mil 

\ 7rm J J pl 

x exp | 

cs * 

1 


, 3 kT G 1 



(61) 


where q c is the condensation coefficient which Young 
argued to be unity and 


Q = 


2(71 — 1 ) L /J 


71 + I R\ Tq \ RiT g 


(62) 


is the nonisothermal correction factor given by 
Kantrowitz in reference 6. 

Transport Properties 

The need to consider transport processes is lim- 
ited to the flow region in which droplet growth 
can occur. The temperature in this region is rela- 
tively low so that viscosity and thermal conductiv- 
ity data are required over a rather narrow tempera- 
ture range. Furthermore, only seven chemical species 
(H 2 O, C0 2 , CO, O 2 , H 2 , N 2 , and NO) are present 
in significant quantities at these temperatures. 

The viscosity of the gaseous mixture can be ob- 
tained by starting with an expression for the viscosity 
of each species in the pure state and then applying a 
mixing rule. The viscosity of each pure species can 
be represented by a polynomial 


Pk — k^G 

n = 1 


(63) 


where the constant coefficients for each species D n ^ 
are given in appendix B. The viscosity of the gaseous 
mixture containing these species can then be ob- 
tained from the expression (ref. 7) 


10 


Pk 


% 10 

fc=1 1 + yr J2 Ylfik/ 
k e=i 


(64) 


where 


4>k,t 


[1 + {Pk/pe) 1/2 i w e/ w k) 1/4 
{8 [1 + m/We)}} 1 / 2 


1 2 


(65) 


Afc = J2 E n,kT£~ l (66) 

n— 1 

where the constant coefficients for each species E n j c 
are given in appendix B. The thermal conductivity of 
the gaseous mixture can then be obtained by using 
the expression (ref. 8) 


10 


1 


Wj2Y t X h+ - „ 

1=1 W E Yt/h 


k= 1 


(67) 


Finally, the equations for the mean free path, 
Prandtl number, and Knudsen number are 


P 

(68) 

p r _ C P^ 

A 

(69) 

Kn = i- 
2r 

(70) 

Kn* = — 
2 r* 

(71) 


and 


where the specific heat at constant pressure for water 
vapor is used to define the Prandtl number. 

Droplet Growth Rate 

The integrated droplet growth equation given by 
Young for droplets with z = — > 1.1 can be written 
as 


6 3 

0 + 1 


lnS±i±? + 


Zj + 6 


+ +!1 ) 


In 


fj+1 ~ 1 

Zj -1 


+ (n + 1 — 9)(zj + i — Zj) 
+ \ { z j+ 1 - z j) = AAt 


(72) 


^ The constant of 1.5 used by Young is somewhat smaller 
than 3 yj 7r/8 = 1.88 from simplified kinetic theory and some* 

nn 8 


what larger than g 4 ^ Q — 1.26 from a considerably more elab- 
orate calculation. 
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This equation also applies in the case with a carrier 
gas present when fi is defined as 


f(l - i^Kn* 

p77 


(73) 


with a carrier gas present when ip is defined as 


= 


fp 


71 + 1 c P) i 1 


Tj ^2 ttRT g 2 ~n. l pl i 


■[T a (pi)-T G ) (80) 


where 


f = y 1 1 — ) 

\ m / 


/m c \ 1/2 

(7c + 1)71 C p,c 

V m / 

(71 + l)7cC p ,i 


The derivation of the factor / is given in appendix A. 

Expressions for the other parameters 9, A, and 
v are 


8 = 2/3Kn* 


A = 


At s (pi)-t g } 

LpltI 


\/ 8 7r 271 
1.5 7i + l 


and 


(75) 

(76) 

(77) 


i/ = 


^1^3 (Pi) 

L 



2 - q c 7i + 1 c p+ t s(pi) 
2 <? c 271 T 


x 



(78) 


The saturation temperature T$(pi) is based on the 
partial pressure of water vapor in the mixture, and 
71 and c P) i are for water vapor. Also note that the 
expression for v contains a term + 1/2 rather than 
— 1/2 given by Young. This difference is believed to 
be due to a typographical error in Young’s paper. 

The integrated droplet growth equation given by 
Young for droplets with z = ^ < 1.1 is 


Zj+ 1 = 1 + T *^ -(zj - 1 ) exp(^Ai) 
r *d+l 

1 ff 7* 

~ 777 exp(t/>AJ - 1) (79) 


where / and v are given by equations (74) and (78). 

The temperature of the liquid droplets is also re- 
quired. Young has developed the following equation 
for Ti which depends on Tq, T s (pi), r*, and r: 

T rp (1 -r*/r)[T s (pi) -Tq) 

T L = T G + (81) 

where 


C £Kn/Pr 

” / 1 < Kn 

1 + 2/3Kn ^ PT 

and £ is given by equation (77). It should be 
noted that equations (72) and (79) imply that the 
droplets retain their identity as they grow and do 
not agglomerate. 

Entropy Production Equation 

The approach taken by Young was to replace the 
differential form of the momentum equation with an 
expression for the increase in entropy due to wa- 
ter condensation. The differential change in en- 
tropy given by Young due to liquid water formation 
also applies to a mixture containing a carrier gas if 
expressed as 



ds — [L 


c p,l Ps(Pl) - T G ]} 



1 

TAP 1 ). 


dw 


(83) 

where the saturation temperature T s (pi) is based on 
the partial pressure of water vapor in the mixture p \ , 
the specific heat of water vapor c p i is used, and the 
latent heat L is evaluated at T 5 (pi). 


Numerical Solution 

The foregoing set of equations can now be used to 
obtain a numerical solution for the adiabatic expan- 
sion of combustion products with nucleation and wa- 
ter droplet growth. The numerical solution is divided 
into three successive parts as indicated in figure 1 . 


Part 1 of Numerical Solution — Stagnation 
Conditions and Mass Flow Rate 


where the subscripts j and j + 1 denote adjacent loca- 
tions in the flow direction and A t is the correspond- 
ing time interval. This equation applies to the case 


Part 1 of the numerical solution involves the spec- 
ification of the feed conditions to the combustor, the 
determination of the adiabatic flame temperature, 
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the computation of the chemical composition of com- 
bustion products in the combustor at a stagnation 
temperature that allows for heat loss, and the deter- 
mination of the mass flow rate. 


Y 2 = a 2,C Y C 

(86) 

£ 

II 

£ 

1 

(87) 


Feed conditions . The first step in computing the 
stagnation conditions in the combustion chamber is 
to define the hydrogen- to-car bon ratio ??h,C a *id the 
number of atoms of carbon in a molecule of fuel n; 
the nitrogen-to-oxygen ratio ^N,Oi the fuel equiva- 
lence ratio <£; the feed temperatures of the fuel, air, 
and oxygen to the combustion chamber, Tjr, T a i r , 
and 7 q 2 ; and the pressure in the combustion cham- 
ber, p Q . The selection of the fuel and the feed tem- 
peratures determines the specific enthalpy of the fuel 
hp as well as and h^ 2 . 

The feed composition in terms of Yp , Yq 2 > and 
>N 2 may be found by first solving equation (16) for 
*o 2 and then equations (13) and (14) for Yp and 

Yn 2 - The specific total enthalpy of the feed can then 
be computed from the equation 


*4 


Yp - Yc - Yj — y 2 
2 


( 88 ) 



(89) 


Y = 


Yu + Yq + y N + Yq — Y\ — r 2 


(90) 


If equations (25) through (30) are now also used 
along with equations (31) through (34), the following 
set of equations can be used to determine the chem- 
ical composition for all 10 species through a series of 
successive approximations: 


y 3 = y 2 


//r P , 2 Tp°\ 1/2 

In p) 


(91) 


ho = Yphp + Yq 2 Hq 2 + Fn 2 £ N2 (84) 

Adiabatic flame temperature . The adiabatic flame 
temperature and the corresponding chemical equi- 
librium composition are determined by the simulta- 
neous solution of the elemental balance expressions 
(eqs. (31) through (34)), the equilibrium relations 
(eqs. (25) through (30)), and the specific enthalpy 
of the mixture with Hq = h Q (eq. (42)). Auxiliary 
equations for the equilibrium constants (eq. (35)) and 
the elemental constants (eqs. (8) through (11)) are 
also required. A solution is obtained by iteration for 
flame temperature Tj . 

Equations (25) through (34) are solved simul- 
taneously by the method of successive approxima- 
tions. For flame temperatures up to 2000 K and 
fuel lean conditions, the major chemical species are 
H 2 O, CO 2 , CO, O 2 , and N 2 . Although significant 
amounts of the other five species may be present, 
estimates for the five major species noted can be 
used for a reasonable initial approximation. If an 
initial estimate of the ratio of the number of mols of 
CO 2 to total mols of carbon in the mixture c^c is 
made, and equations (31) through (34) are used with 
Y$ = Yj = Yg = Yq = Y\q = 0, it follows that a first 
approximation for the mol numbers of the five major 
species and the total mol number of the mixture is 



(85) 


Ts = 


Y\Y 3 


K p ,iY 2 W 

Yio= (Kp^Ye ) 1 ' 2 (93) 

r 6 = 10 (94) 

Yg = (K p> 3 Y 4 Y 5 ) 1/2 (95) 

Y 7 = \K Pt 4 Y 5 Y^j (96) 

T 8 =(ifp )5 T 4 Y^) 1/2 (97) 

Yi = Yh-2Y s -Y 7 -Yq {w) 

Y 2 = Y c - Y 3 (99) 

y 4 = ^-y c -yi-y,-y,-y,-yio (100) 


Y 4 


[y 4 + Yo-Yc-n-Y2-Y 8 -Y9-Y 10 


10 

y = £n 

k=l 


(101) 

( 102 ) 


It should be noted that equation (100) led to step- 
to-step oscillations in the solution for higher tem- 
peratures. In order to avoid this problem, equa- 
tion (100) was rewritten as equation (101) to pro- 
vide some damping. The numerical value of Y 4 com- 
puted in the previous iteration step is used in the 
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right-hand side of equation (101). After a stable so- 
lution for the mol numbers Y ^ for k = 1 to 10 is 
achieved for a given temperature, the correspond- 
ing specific enthalpy of the mixture Hq is computed 
from equation (42) and compared with the known 
total feed enthalpy h Q . If Hq turns out to be less 
than h Qy the next approximation to the temperature 
is increased and if greater, decreased. The final solu- 
tion for the adiabatic flame temperature is obtained 
when Hq — h Q . 

Stagnation conditions . This study does not at- 
tempt to compute the heat loss from the combustor, 
however, the influence of an arbitrary heat loss has 
been taken into account by assuming a stagnation 
temperature T 0 that is somewhat less than the adia- 
batic flame temperature. The chemical composition 
at T 0 can then be determined directly by using the 
scheme just presented. The stagnation enthalpy h Q 
and the stagnation entropy s Q are then computed 
from equations (42) and (44). 

Mass flow rate. The mass flow rate is determined 
in the following way. An isentropic expansion from 
the combustion chamber is assumed and the mass 
flux pqU computed along the nozzle to find the maxi- 
mum value of pqU. Since h 0 and s are constant along 
the expansion, all flow properties can be computed 
at a selected temperature that is less than T 0 . This 
computation uses an iteration on the pressure. 

The choice of an approximate pressure, say p', at 
the selected temperature T allows the computation 
of a chemical equilibrium composition. This in turn 
permits the direct computation of the density by 
equations (45) and (46). The resulting density pq 
and composition can then be used in equation (47) 
to compute a pressure p. A comparison between p f 
and p indicates the need to continue the iteration 
process. When the iteration scheme gives a value of 
p sufficiently close to p f , the pressure is determined. 
The flow velocity U at this state is then computed 
from equation (3) with no liquid water present; that 
is, 


U = [2(ho - h G )} 1/2 (103) 

The mass flux pqU can therefore be found at 
selected temperatures corresponding to successive 
locations along the nozzle. The maximum mass flux 
PqU* occurs at the nozzle throat where the nozzle 
cross-sectional area is A * . The mass flow rate is given 
then by equation (1) with no liquid water present; 
that is, 


m = p* G U*A* (104) 


This is the mass flow rate of the total mixture along 
the entire length of the nozzle including the down- 
stream conditions with liquid water formation. 

Part 2 of Numerical Solution — Isentropic 
Expansion and Beginning of Nucleation 

Part 2 of the numerical solution involves the com- 
putation from the nozzle throat to the beginning of 
nucleation. It is assumed that the flow expands isen- 
tropically and remains in chemical equilibrium with 
no liquid water formation. At each step along the 
nozzle the nucleation rate is computed to determine 
at what point to include nucleation and liquid droplet 
growth. 

Isentropic expansion . The isentropic expansion 
from the nozzle throat is obtained by noting that 
s = s Q and h Q = Constant along the nozzle. Start- 
ing at T* and p*, a temperature T less than T* is 
selected. The pressure is then determined by itera- 
tion and the values of p< 7 , hQ , and U computed by 
the scheme given in part 1 of the numerical solution. 
The area of the nozzle can then be computed from 
equation ( 1 ) with w = 0 ; that is, 

P G U 

The location of this particular point along the nozzle 
is then determined by using the relationship that 
defines the cross-sectional area as a function of nozzle 
position x. This computation in general requires the 
inversion of the expression of A as a function of x or 
an iteration to find x for a given value of A. 

Beginning of nucleation * The nucleation rate J is 
computed from equation (61) where the appropriate 
values of 71 , Hq j, L, p 0 c , <r, and r* are obtained 
from equations (39), (43), (56), (57), (58), and (59), 
respectively. This computed value of J is compared 
with a value J m \ n chosen such that a smaller value of 
J does not generate enough nuclei to affect the flow 
throughout the nozzle. Successively lower values of 
T are taken along the nozzle until J > J m \ n < At 
this point and downstream, the number of nuclei 
formed is taken into account, and the droplet growth 
equations are used. 

Part 3 of Numerical Solution — Nucleation and 
Droplet Growth 

Part 3 of the numerical solution deals with the 
computation of the nucleation process to determine 
the number of droplets formed at successive loca- 
tions along the nozzle and the subsequent growth 
of the liquid water droplets. The starting point is 
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the condition determined in part 2 of the numeri- 
cal solution where the nucleation rate just exceeds 
a threshold value. The index j is used to denote a 
location along the nozzle, starting with j — 1 at the 
beginning of part 3 (fig. 1), where all properties are 
known. This also includes the chemical composition 
of the gas phase. Since essentially no further chemi- 
cal reaction takes place once nucleation and droplet 
growth occur, the only change in the gaseous mix- 
ture composition is due to the formation of liquid 
water. The appropriate expressions for Yj, Y*. (with 
k = 2, 3, ..., 6), and W at any position along 
the nozzle are thus given by equations (17), (18), 
and (19). The computation proceeds in general at 
location j where T , p, p, t/i, f/, s, A, z, and w are 
known. All other required parameters at j can be 
computed from this set. Note also that all properties 
at j — 1 are known. 

Computational steps along nozzle and nozzle area . 

The first step in the determination of properties at 
j + 1 is to increase xj by Ax so that 

Xj+i=Xj + Ax (106) 

and to determine the nozzle area from the expression 

Aj + 1 = f{x i+ 1 ) (107) 

An iteration scheme based on adjusting Tj+\ as the 
primary variable is then used to find a solution for 
all properties at j -j- 1. Secondary iterations within 
the loop for Tj + i are also required. 


Although the droplet size of this set formed at j will 
change as it moves downstream, the number An 1 
in the set will remain constant unless some of the 
droplets evaporate. Agglomeration of droplets is not 
taken into account in these calculations. Note that a 
set of droplets characterized by a single index i will 
be formed at successive positions along the nozzle. 
At any position j, there could be as many as j sets 
of droplets each with a radius of rj- in the amount of 

An 1 droplets per mass of total mixture. 


Droplet growth . The computation of the droplet 
radius at j + 1 employs equation (72) and appropri- 

ate auxiliary equations for z\ = -r 2 - > 1.1 or equa- 

3 ' * ,j 1 

r « 

tion (79) and auxiliary equations for z l ■ = < 1 . 1 . 

3 r +J 

Either equation requires the time interval that cor- 
responds to the step size Ax; that is, 


It should be noted that the derivation of equa- 
tion (72) by Young used the argument that the prop- 
erties of the gaseous mixture do not change between 
locations j and j+ 1. This means that equation (72) 
can be used to compute rj- +1 with only knowledge 
of properties at j . The solution for droplets in the 
size range of z x >1.1 can be obtained by rewriting 
equation (72) as 


Gas temperature . The first approximation for the 
gas temperature at j + 1 is T J+1 = Tj and it is further 
assumed that initially Pj+\ = Pj and Uj+ \ = Uj. 
The value of Tj+i is adjusted as a result of each 
iteration step, and Tj +1 = Tj +l is defined at this 
point and used to test for convergence. 

Formation of new droplets . The nucleation rate is 
assumed to vary as exp(fcx) as suggested by Young, 
so that the appropriate mean nucleation rate ex- 
pressed in terms of new droplets formed per unit vol- 
ume per unit time at j over the volume from j — i 
to j -f j is the logarithmic mean value; that is, 



J j± 1/2 ~ J j- 


1/2 


(108) 


where J an d J j_\j 2 are computed from equa- 
tion (61). The number of new droplets formed at j 
per unit mass of total mixture is therefore 


An l “ J 


JjAjAx 

m 


(109) 


F = Zj +1 - Zj - AAt (111) 

where 


+ («ir+ n )">(^-i) 

+ (n + i -e)z) + ^ z f (ii2) 

The application of a Newton iteration scheme re- 
quires dF/dz l j +1 which is 

dF _ e 3 

" (» + i) (*;■+,+*) 

I 1 + 0(0 + 1 ) 

+ (n + i ~ 
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(113) 



so that the nth iteration gives 


F(n) 


z l j+1 (n + l)-z l j+l (n) dF{n)/dzl 


(114) 


i+i 


The droplet radius at j + 1 when the iteration has 
been completed is 


r l 

1 


= z 




(115) 


Equation (79) is used for the smaller droplets 
when < 1.1- In this case, the solution for z hi 
requires knowledge of 1 as well as properties at 
j. An expression for dr*/dt is 


±1 = ~ '*•> (116) 
dt At v ; 

It should be noted that Young’s analysis assumed 
that tp in equation (79) and given by equation (80) 
is a constant over the interval from j to j + 1 but 
that the solution for from equation (79) requires 
knowledge of r+j+i. This means that the computa- 
tion of rj. +1 when z*- < 1.1 must use successively 
revised values of r+j+i that correspond to the suc- 
cessively revised values of Tj+\. 


Mass of liquid water . The mass of liquid water 
droplets that have a radius r j+\ is g7r pi 
so that the mass of liquid water per total mass of 
mixture for all droplets of type i at j + 1 is 

Wj + 1 = !*7>L ( r }+l) 3 An * ( 117 ) 

The total mass of liquid water per total mass of 
mixture at j + 1 is 


i+1 

w j+\ = w j+i ( 118 ) 

i=l 

Enthalpy and entropy of liquid water. Now in order 
to compute the enthalpy and entropy of the liquid 
water in the mixture, the temperature of the liquid 
droplets of each type i , must be determined 

from equation (81). The enthalpy of the liquid at 
y + 1 can then be computed from equation (50) using 
T l L - + i to get h { L j+l and h L j+i is obtained from 


The entropy of the liquid water is computed like- 
wise from equation (54) using * +1 to get s'lj+i 
so that 


•w+i = e«$ + i*w+i (‘20) 

i 

Entropy production. The next step is to compute 
the entropy increase due to liquid water formation 
by equation (83) with ds replaced by Asj and dw by 
A Wj where 

Awj — Wj + 1 — wj (121) 

The temperature and pressure used to determine 


properties in equation (83) are 


t gj+i + t gj 
Tg= 2 

(122) 

and 


Pj + 1 + Pj 
P ~ 2 

(123) 


The value of L is evaluated at T s (pi), and C P,1 is 
a very weak function of temperature so it can be 
evaluated at Tq. The entropy of the total mixture 
can therefore be written as 

sy+1 = Sj + Asj (124) 

Two-level iteration to determine T^. At this point, 
the iteration scheme proceeds with an attempt to 
find a value of Tqj+i that satisfies the independent 
computation of sqj+ i- Simultaneously an iteration 
is carried out to determine a value of PGj+1 by using 
the known area at j + 1, Aj+ 1 , and the mass flow 
equation. The two iteration schemes are as follows: 

1. Iteration using entropy of gas mixture: One 
equation for the entropy of the gaseous mixture can 
be written as 


n _ 1 ~ s Lj + l 

GJ+l i _ Wj . +1 

where all terms on the right-hand side have been 
computed. A second and independent equation is 
given by equation (44) from which S GJ + 1 is com - 
puted using PGJ + 1 an( l ^GJ- hi* The iteration at this 
point assumes that s qj+ l correc t> fixes PGJ+U 
and adjusts Tqj+ \ until 

(4 ~ fc)j±l < e (126) 

Cp R 

where 

k 
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If this condition is not satisfied, the value of Tqj+ i 
is revised according to the expression 


T j+ i( n + 1) = i( n ) 


( I * * * S G 


S G )i + 1 


Cp - R 


(127) 


where T J+1 = T Gj+1 . 

2. Iteration using nozzle area: After the con- 
dition in equation (126) is satisfied, Ty +1 is used in 
equation (42) to compute h f G ^^. The specific en- 
thalpy of the total mixture is then computed from 


h 'j+ 1 — C 1 _ w 3 + l)h' G j + l + ^Lj'+l (128) 


The value of U at j + 1 is then computed from 


U J+1 = [2(h 0 - h' j+1 )] U2 (129) 

A value for the cross-sectional area of the nozzle at 
j H- 1 is then computed from equation (1); that is, 


, m ( 1 — Wj+ 1 ) 

>+1 PGJ + lU j+ i 

and compared with the known area Aj+\ by using 
the expression 



^j‘+l — A? + 1 
Aj + 1 

If this condition is not satisfied, the density of the 
gaseous mixture is computed from the equation 



PGj+1 = 


m(l -wj+i) 

Aj+iUj+i 


(132) 


and the iteration process returns to the computation 
°f s gj+i usm E ^is value of Pcj+1 and the last 
computed value of Tqj+i. The iteration continues 
until the condition in equation (131) is satisfied after 
which Pj+i is computed from equation (47). 

The iteration for now includes the test that 


\Tj+l 


I 

x j + ll 


T? 

1 j + 1 


< £ 


(133) 


If this test is not satisfied, the iteration process 
returns to the beginning (gas temperature) with T J + i 
equal to its last value. When this test is satisfied j is 
increased by 1 with x increased by Ax (step along the 
nozzle and nozzle area). This process is continued to 
the nozzle exit. 


Computer Program 

A computer program, FIRACON, has been de- 
veloped to solve the finite-rate condensation problem 
using the solution technique described in the present 
paper. The program is written in FORTRAN Ver- 
sion 5 language for the CDC® CYBER series digital 
computer system, model 860, with Network Oper- 
ating System. The program requires approximately 
134000 octal locations of core storage, and a typi- 
cal case requires approximately 200 seconds of cen- 
tral processing unit time. Program input require- 
ments and operating instructions are summarized in 
the code which is included as appendix C. 

Application of Computation Scheme to 
Langley 8-Foot High-Temperature Tunnel 

The numerical scheme presented in the previous 
section has been applied to the Langley 8-Foot High- 
Temperature Tunnel (8'HTT) for a selected set of op- 
erating conditions. The 8'HTT is a large combustion- 
driven wind tunnel that burns natural gas (mostly 
methane) and air at high pressure in a cylindrical 
combustor approximately 1 m in diameter and 7 m 
long. The combustion products which are primar- 
ily H 2 O, C0 2 , O 2 , and N 2 with smaller amounts of 
H 2 , CO, H, O, NO, and OH, expand as the test 
gas from the combustor through an axisymmetric 
conical-contoured nozzle to a test section with a di- 
ameter of approximately 2.4 m (8 ft). The nozzle 
length is approximately 16 m. 

Numerical results are presented for five cases that 
correspond to operating conditions of the 8'HTT. 
The combustor temperature, pressure, and fuel 
equivalence ratio for these cases are given in table 1. 
The mole fraction of water vapor formed in the com- 
bustor for each case is also listed. The conditions 
for case 1 were established by selecting a combustor 
pressure of 50 bars and an equivalence ratio that gave 
a flame temperature of 2000 K. An arbitrary heat 
loss from the combustor was included by selecting a 
stagnation temperature that was 100 K less than the 
flame temperature, that is, 1900 K. Conditions for 
case 2 were obtained likewise but for a pressure of 
250 bars. Conditions for case 3 and case 4 were ob- 
tained in a similar fashion but for a temperature of 
1600 K. Conditions for case 5 were obtained for T 0 
— 1900 K and p Q = 250 bars, the same as for case 2, 
but for additional oxygen and an equivalence ratio 
that yielded a mole fraction of 0.21 for O 2 in the test 
stream. Note that the mole fraction of water vapor 
in the combustor is 0.154 for cases 1 and 2, 0.122 for 
cases 3 and 4, and 0.156 for case 5. The values used 
for the coefficients g c , a, a c , /?, and J m \ n are also 
indicated in table 1. 
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The nozzle contour of the 8'HTT was used to de- 
fine the variation of the cross-sectional area along 
the nozzle. Tabulated values of selected wall coordi- 
nates r w and computed boundary-layer-displacement 
thicknesses 6* are tabulated in table 2. The numeri- 
cal values of 6* were computed from a modified ver- 
sion of the boundary-layer scheme presented in refer- 
ence 9. The inviscid nozzle contour that was used for 
the quasi-one-dimensional computation was obtained 
by subtracting 6* from r w at each location of x listed 
in table 2. A spline fitting scheme for intermediate 
locations was applied. 

The computed results for case 1 with T 0 — 
1900 K, p 0 = 50 bars, and <f> = 0.798 are shown 
in figure 2. Figure 2(a) shows the difference be- 
tween saturation temperature and the gas temper- 
ature Ts(pi) - Tq ^ the nucleation rate J, and the 
mass fraction of liquid water formed. The tempera- 
ture difference T s (pi) - Tq is a measure of the de- 
gree of supercooling and provides the driving force 
for nucleation and droplet growth. The test gas 
first becomes saturated with respect to water vapor 
when T s {p\) = Tq. This occurs at a nozzle loca- 
tion just beyond x = 6.5 m. The temperature dif- 
ference T s {pi) — Tq then increases to a maximum 
value of about 45 K at approximately x = 12 m after 
which T s {pi)-Tq decreases rapidly. Figure 2(a) also 
shows that the nucleation rate J rises exponentially 
and then decreases rapidly as T s {pi) - Tq decreases. 
The mass fraction of liquid water begins to be signif- 
icant at about the point of the maximum nucleation 
rate. The liquid water then continues to increase to 
the end of the nozzle. 

The formation of liquid water under these condi- 
tions results in entropy production as shown in fig- 
ure 2(b). The flow is isentropic until liquid water is 
formed. The rather rapid formation of liquid water 
initially causes a correspondingly rapid increase in 
entropy followed by a more gradual increase as the 
rate of liquid water formation decreases. This trend 
is to be expected from equation (83). 

The entropy production due to liquid water for- 
mation affects the static pressure distribution along 
the nozzle as shown in figure 2(c). The dashed curve 
indicates the pressure distribution for an isentropic 
expansion if no liquid water was formed. The solid 
curve represents the computation that takes into ac- 
count nucleation and droplet growth. Note that the 
process follows the isentropic (dashed) curve to just 
beyond a nozzle location of x =■ 12 m. At this point 
the static pressure decreases less rapidly than indi- 
cated by the isentropic process. The pressure then 
increases somewhat to a maximum value after which 
it decreases in parallel to the isentropic curve. This 
departure from the isentropic pressure distribution is 


due to the transfer of heat to the gaseous phase from 
the liquid droplets. At the nozzle exit the static pres- 
sure for case 1 with condensation is approximately 
25 percent larger than that computed for isentropic 
flow without condensation. 

The growth history of the liquid droplets is shown 
in figure 2(d). The single curve labeled r* is the 
critical droplet radius computed from equation (59). 
Droplets that are larger than r* grow, whereas those 
less than r* evaporate. The series of curves indi- 
cates the growth of the sets of droplets that were ini- 
tially formed at successive stations along the nozzle. 
Droplet sets that are formed but eventually evapo- 
rate are also shown. Only a portion of the numerous 
sets of droplets formed is included in figure 2(d). 

The computed results for case 2 are given in 
figure 3 and are also for T 0 = 1900 K but at a higher 
pressure of p Q — 250 bars. Figure 3(a) shows that 
T s {pi) “ Tq increases earlier and more rapidly for 
p 0 = 250 bars than shown in figure 2(a) for p 0 = 
50 bars. Case 2 shows a slightly lower maximum 
value of T s {pi) - Tq and also an earlier and more 
rapid decrease in T s (p\) — Tq along the nozzle. The 
nucleation process also begins further upstream for 
the higher pressure case and displays a more rapid 
increase and decay. The maximum nucleation rate 
also achieves a higher value for the higher pressure 
case, Note also that liquid water forms earlier and 
more rapidly for p 0 = 250 bars compared with p Q — 
50 bars. The total amount of liquid water formed for 
the higher pressure case is approximately 60 percent 
more than for case 1. 

Figure 3(b) shows the entropy production for 
case 2. A comparison with figure 2(b) shows that 
the entropy rises more rapidly for the higher pressure 
case. This is due to the more rapid formation of 
liquid water. The entropy also attains a nearly 
level value, but the total entropy increase is less 
for p 0 — 250 bars than for p 0 = 50 bars. This 
smaller entropy increase is a direct result of the 
earlier collapse of the supercooling. 

The static pressure distribution along the nozzle 
for case 2 is shown in figure 3(c) . A comparison of fig- 
ure 3(c) with figure 2(c) shows that the higher pres- 
sure case results in an earlier and more pronounced 
departure from an isentropic process. The static 
pressure at the nozzle exit for case 2 with conden- 
sation is approximately 33 percent larger than that 
computed for isentropic flow without condensation. 

The droplet growth history for case 2 is shown in 
figure 3(d). The results are similar to case 1 except 
the liquid droplets form earlier, grow more rapidly, 
and attain larger sizes. 

The results from case 3 are included to in- 
dicate the effect of a lower temperature and are 
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presented in figure 4. A comparison of these re- 
sults for T 0 = 1600 K with those for case 1 with 
T 0 = 1900 K shows that the lower temperature case 
results in an earlier onset of nucleation and conden- 
sation. More liquid water is formed even though the 
total water vapor initially present in the test gas is 
less for case 3 than for case 1. The static pressure 
distribution plotted in figure 4(c) for case 3 shows a 
more pronounced departure from an isentropic pro- 
cess than was noted for either case 1 or case 2. 
The static pressure at the nozzle exit computed for 
flow with condensation is approximately 57 percent 
greater than that computed for isentropic flow with- 
out condensation as shown in figure 4(c). The re- 
sults for case 4 with T 0 = 1600 K, p Q = 250 bars, 
and <f> = 0.62 are presented in figure 5 and can be 
compared with the other cases. The results for case 5, 
which are for additional oxygen to give a test stream 
containing 0.21 mole fraction of O 2 , are presented in 
figure 6. A comparison of these results with case 2 
(fig. 3), which is also for 1900 K and 250 bars, shows 
very little difference. The somewhat different com- 
position of the combustion products produces an ex- 
pansion that gives a water vapor saturation point for 
case 5 that is slightly downstream of that for case 2. 
All the curves shown in figure 6 for case 5 are there- 
fore shifted slightly downstream relative to case 2. 

The purpose of presenting these five cases is to 
show a set of results that indicates the magnitude 
and the general trends for the effects of water con- 
densation in the Langley 8'HTT. The numerical val- 
ues used for the coefficients q c =1, a = 8, and (3 = 2 
are those suggested by Young for pure steam. The 
value of a c = 1 was used to account for the non- 
condensable carrier gas. Computations with values 
of a c < 1 were also carried out but are not included, 
and they give similar but less pronounced results. 
The value of J m \ n used was adjusted until a smaller 


value had no effect. More appropriate values of the 
empirical coefficients a,a c , and ) 3 might be deter- 
mined from experimental data. These data would 
include static pressure distributions along the nozzle 
and water droplet size and concentration measure- 
ments for supersonic expanding flow of combustion 
products. 

Concluding Remarks 

An analysis and numerical scheme has been de- 
veloped to treat the supersonic expansion of com- 
bustion products that takes into account nucleation 
and finite-rate growth of liquid water droplets. The 
scheme has two limitations: the flow is assumed to 
be quasi-one-dimensional and empirical coefficients 
are required. This analysis has specific application 
to the computation of flow properties in combustion- 
heated wind tunnels such as the Langley 8-Foot High- 
Temperature Tunnel (8'HTT). 

Sample computations included in this paper are 
based on the nozzle contour of the Langley S^HTT. 
The numerical values used for the empirical coeffi- 
cients in these calculations are those suggested by 
Young for pure steam. The results indicate that 
the free-stream static pressure can be significantly 
greater than that computed for isentropic nozzle flow 
without condensation. The computed entropy pro- 
duction also indicates a loss in total pressure relative 
to isentropic flow. These computed results suggest 
that this scheme can be used as a tool to interpret the 
calibration and flow measurements in combustion- 
heated wind tunnels. A calibration of combustion- 
heated wind tunnels must take the effects of water 
condensation into account. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
July 29, 1988 



Appendix A 

Extension of Young’s Droplet Growth 
Equations To Include Carrier Gas 

This appendix extends the droplet growth equa- 
tions given by Young in references 3 and 4 for pure 
steam to the case of water vapor in the presence of 
a carrier gas. The carrier gas takes into account 
the noncondensable gases contained in the products 
of combustion. The equations for mass and heat 
transfer given by Young for water vapor are used 
and the appropriate equations for a carrier gas are 
introduced. 

Young’s analysis begins by writing the equations 
for the free molecular transport regime and for the 
continuum transport regime and then combines these 
equations to find droplet growth equations that ap- 
ply to the intermediate regime. Young postulates 
that the transport processes can be described by free 
molecular transport from the droplet surface out to a 
distance of a few mean free paths and by continuum 
transport beyond. The interface between these two 
regimes is taken to be at r + (3t where (3 is a parame- 
ter of the order of unity. Young has chosen a value of 
1 3 = 2 . At this interface the pressure and tempera- 
ture are denoted by pi and T±. 

Free Molecular Transport Regime (Kn > 1) 

The mass transfer equation for water in the free 
molecular regime from the droplet surface out to the 
interface given by Young is 


where p c { is the partial pressure of the carrier gas at 
the interface, R c = RjW Cy and a c is the thermal ac- 
commodation coefficient for the interaction between 
the carrier gas and a water droplet. It should be 
noted that the thermal accommodation coefficient for 
water vapor with a water droplet was assumed to be 
unity by Young as is implied by equation (A2). The 
total energy transfer in the free molecular regime due 
to both water vapor and carrier gas is then the sum 
of equations (A2) and (A3); that is, 


L 


dM 

dt 


= 47rr 2 


Pi 71 + 1 

\J 2iiRTi 271 


c P ,i(T L -Ti)f 


(A4) 


where 


/ = yi 




{lc + 1)71 Cp,c 
(71 + lbcC p ,i 


OL C 


(A5) 

This expression for / is equation (74). Note that / 
reduces to unity for the pure steam case. 


Continuum Transport Regime (Kn < 1) 


Young points out that even when the conden- 
sation rate is high, the pressure drop beyond the 
interface due to bulk flow of water vapor towards 
the droplet is negligible so that 


Pi = P 


(A6) 


dM 2 q c 2 

— — = 47rr 

dt 2 — q c 


Phi 


<?e Ps{TL)T) 


m yj2ir R\T{ q ° \j2^R\Tj Jm 


(Al) 


and the energy transfer equation in this regime is 


The energy transfer equation in the continuum 
regime given by Young is 

L -jr — 4 *{r + /3£) 2 -A-,(T l - 7b) (A7) 

at r + 0t 



= 4nr 2 


Pl,« 

JznTUTi 


71-kl 

271 


c pA( t L “ T i) 


(A2) 

The mass transfer equation is appropriate as it stands 
with pi j representing the partial pressure of water 


vapor at the interface and R\ = R/W\. Equa- 
tion (A2), however, includes only the energy transfer 
due to water. An additional equation is required to 
account for the parallel energy transfer due to the 
carrier gas. The equation for the energy transfer due 
to the carrier gas has the same form as equation (A2); 
that is, 



= 47r r 2 


Pc, i 7c + 1 

faRcTi 2lc 


(A3) 


Equation (A6) can be used as it stands and so 
can equation (A7) if A is taken to be the thermal 
conductivity of the gas mixture. 

Intermediate Regime 

The mass transfer equation in the intermediate 
regime is obtained by combining equations (Al) and 
(A6); that is, 


dM 

dt 


2 q c 

2 Qc 


47T r 2 



Qe Ps{T L) r) 
Qc \/2irRiT L 


(A8) 

The energy transfer equation in the intermediate 
regime can be developed by replacing pi in equa- 
tion (A4) with p and combining that result with 
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equation (A7) to obtain 



r 


47rr 2 (A/r)(T L -r G ) 

1 | ^ 2~n 1 

r + /3t r^p Tl + 1 c p ,l/ 


(A9) 


Now if equation (68) is used to eliminate p and 
equations (69) and (70) are used for Pr and Kn, it 
follows that 



47rr 2 (A/r)(r L -T G ) 

1 £(Kn/Pr)V7y7b 

1 + 2/9Kn + f 


(A 10a) 


and with the additional approximation that \/T x /Tq 
can be taken as unity, 

dM _ 47rr 2 (A/r)(r i -T G ) 

* "[r+W^j 
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Young continues his analysis by developing a 
relationship between T^—Tq and T s (pi) — Tq. First, 


{T L - T G ) = [T.Cpi) - t g ] - (t' s - T L ) 

+ [Ti - T a ( Pl )} (A12) 

where T f s is the saturation temperature that corre- 
sponds to p 5 (T^,r). Young states that a very good 
approximation to the second term is 


i T 's - T L ) = r -y[T s {pi) - T g ] (A13) 

The third term, T f s — T s (pi), is obtained in the 
following way. Equation (A4) with p{ = p and 
equation (A8) are used to eliminate dM/dt to obtain 
an expression for p 5 (T^,r)/pi; that is, 


Ps{T L ,r) _ qc_ 
Pi 7e 



2- q c 7l + l Cp 7 1 
2 q c 2^1 L 


An expression for - T x in equation (A 14) can be 
obtained by equating equation (A4) with p x — p to 
equation (A 10a) to get 


where 


Tl - Ti = 6 (T l - T G ) 

fKn/Pr 

JyjTdTi , £Kn 
T + 2/3Kn + Pr 


6 = 


(A16) 


(A 17a) 


T he add itional approximation used by Young that 
VTb/Ti is unity leads to 


£Kn/Pr 


+ 


Kn 


(A17b) 


which is equation (82). Young also postulates that 


9 c _ i a(r L - Tj) 
9e T s (pi) 


(A18) 


where a is a constant and higher order terms are 
neglected. Also note that 

-(^-f <-> 

Substituting equations (A18) and (A19) into equa- 
tion (A14), taking the logarithm of both sides of the 
resulting equation, and then using the approximation 
In (1 + x) — x for x 1 yield 

In Fs(^’ r ) _ q CTl ~ Ti) \{T l - Ti) _ 2 - q c 71 + 1 
Pi T s (pi) Ti 2 q c 271 

>t < A20 > 

Substituting this expression into equation (A 15) and 
using equation (A 16) give 



T' s -T s { P i) = u6{T l -T g ) (A21) 


x {Tl - Ti) 


J_ / m 

S/l Vmi . 


, 1/2 


(A14) 


where 


Young then uses an approximate integral form of the 
Clausius-Clapeyron equation, 


TJ _ RlTsiPi) [ , lT' 5 (pi) 
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The additional approximation that T s (pi)/T{ is unity 
leads to 

RiTsipi) r . 1 2 -g c 71 + I c Py iT s (pi) 

u = L L 2 ~ ~2^T “ 2 ^ I 

X Z(^) 1/2 1 (A22b) 

Vi \mx / J 

which is equation (78). A relationship between 
(Tl - Tq) and [T 5 (pi) - Tq] can now be formed by 
using equations (A 13) and (A21) in equation (A 12) 
to obtain 


[t l - t g ) = [ Ts (pi) " T d (A23) 


This expression for T \ — Tq can be substituted into 
equation (A 10b) and since the mass of a liquid water 
droplet is 

M = ^nr 3 p L (A24) 

the droplet growth equation becomes 


dr 

dt 


A(1 - r m /r)[T a (pi) - Tq] 


(W)[r^ + ^ P f - 


(A25) 


This equation reduces to that given by Young for 
pure steam for / = 1. 

For droplets with z = ~ > 1.1, Young defines a 
Knudsen number based on a critical droplet radius 
as 


Kn. = (A26) 

2r* 

and assumes that the vapor properties remain con- 
stant over an integration step so that equation (A25) 


can be written as 


dz 

dt 


A(i-i/z)[r s (pi)-7b] 


r , 9 [ 2 2 , $(Kn*/Pr)(l - u) 

L PL r i [z + 2/3Kn* + 


(A27) 

This equation can be integrated over a time step At 


from Zj to Zj+i to give equation (72) where the terms 
n, /, 0, A, and v are defined by equations (73) 
through (78). 

For droplets with r/r* < 1.1, Young uses the en- 
ergy transfer equation appropriate in the free molec- 
ular regime; that is, equation (A4) with p t — p and 
T{ = Tq. Again if equation (A23) but with 5 = 1 and 
equation (A24) are used, it follows that for £ < 1.1, 


dt p L L\j2nRT G 


Now define 


(A28) 


f : 


ryj2'KRTc 




so that 


l = r *(i-7) 

*/dt 
d(r - r*) 


(A29) 

« -V n (A30) 

Subtraction of dr * /dt from both sides of this equation 
gives 


dr* 


dt 


= (A31) 


At this point, Young makes the approximation that 
ip is constant over an integration step so that equa- 
tion (A31) can be integrated over a time interval At 
to give equation (79). An expression for ip given by 
equation (80) is obtained from equation (A29) with 
r set equal to ry, and all other quantities in the ex- 
pression for ip are also evaluated at station j . 
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Appendix B 

Thermodynamic and Transport Property 
Data for Products of Combustion 

This appendix lists a set of coefficients from 
which thermodynamic and transport properties can 
be computed for products of combustion. The 
thermodynamic properties can be computed over the 
temperature range from 200 K to 3000 K and the 
transport properties up to 500 K. The coefficients for 
the thermodynamic properties were obtained from a 
least-squares fit to tabulated heat capacity and equi- 
librium constant values given in the JANAF Thermo- 
chemical Tables (ref. 10). The coefficients for com- 
puting viscosity were determined from a least-squares 


fit to tabulated data presented by Touloukian, Sax- 
ena, and Hestermans (ref. 7) and the coefficients 
for computing thermal conductivity using tabulated 
data presented by Touloukian, Liley, and Saxena 
(ref. 8). 

Table B1 lists coefficients A n % (n= 1, 2, .... 6) 
that appear in equation (36) for the heat capac- 
ity, the integration constant for the enthalpy Ay 
that appears in equation (42), and the integration 
constant for the entropy Ag ^ that appears in equa- 
tion (44). The numerical values listed correspond to 
units for C p ^ in calories per mol-kelvin. The coef- 
ficients listed in table Bl(a) apply to a temperature 
range of 200 K < T < 1000 K and table Bl(b) apply 
to 1000 K < T < 3000 K. 


Table Bl. Heat Capacity Coefficients A n ^ (n = 1, 2, ..., 6) and Constants of Integration 
for Enthalpy Ay ^ and Entropy Ag 


k 



A n,k 

For n = 



A 7 ,k 

A 8,k 

1 

2 

3 

4 

5 

6 

(a) 200 K < T < 1000 K 

1 

6.41000E+01 

7.59151E+00 

— 8.41829E— 04 

6.68750E-06 

— 5.19580E— 09 

1.55769E-12 

— 6.04359E+04 

2.06410E+00 

2 

1.04750E+02 

3.37574E+00 

2.38924E-02 

-2.68777E-05 

1.71801E— 08 

— 4.69551E- 12 

-9.65120E+04 

2.61229E+01 

3 

-1.60500E+02 

9.06284E+00 

— 9.91073E-03 

1.97943E-05 

— 1.49575E— 08 

4.10256E- 12 

— 2.79118E+04 

-2.76320E+00 

4 

4.94000E+01 

6.94324E+00 

— 3.59602E — 03 

1.51603E-05 

— 1.54277E-08 

5.20513E— 12 

— 2.29792E+03 

1.01347E+01 

5 

— 4.26900E+02 

9.83533E+00 

— 7.22941E— 03 

8.48001E— 06 

— 4.45280E— 09 

1.01282E— 12 

— 2.44860E+02 

-2.44460E+01 

6 

-1.51550E+02 

8.84079E+00 

-8.29361E-03 

1.52780E— 05 

— 1.03811E— 08 

2.52244E— 12 

-1.51928E+03 

— 3.22920E+00 

7 

0 

4.968 OOE+OO 

0 

0 

0 

0 

5.06219E+04 

— 9.14000E— 01 

8 

1.06200E+02 

5.01442E+00 

— 8.54283E— 04 

1.7432 IE -06 

— 1.46562E— 09 

4.55128E — 13 

5.74789E+04 

1.04430E+01 

9 

1.65300E+02 

6.18389E+00 

2.60375E— 03 

— 5.47821E— 06 

5.88287E— 09 

-2.02564E— 12 

6.45433E+03 

8.62130E+00 

10 

— 9.10512E+00 

8.49215E+00 

— 9.77414E— 03 

2.32087E — 05 

— 1.98449E— 08 

6.05054E— 12 

1.93658E+04 

3.97820E+00 

(b) 1000 K < T < 3000 K 

1 

3.49199E+03 

-4.64184E+00 

1.68818E — 02 

— 7.39530E — 06 

1.68507E— 09 

— 1.58341E— 13 

— 7.74261E+04 

7.74501E+01 

2 

-2.22684E+03 

1.37397E+01 

2.55399E— 03 

— 1.4250 7E — 06 

3.76285E—10 

— 3.87790E— 14 

— 8.53153E+04 

-3.47509E+01 

3 

—3. 5834 IE -|- 02 

6.31771E+00 

3.26443E— 03 

— 1.66622E— 06 

4.13442E— 10 

— 4.03866E— 14 

-2.62483E+04 

9.46980E+00 

4 

— 1.5392 2E+ 03 

1.07380E+01 

— 1.38694E— 03 

5.81156E— 07 

— 5.58070E-11 

— 3.284 50E— 15 

5.83496E+03 

-L64089E+01 

5 

1.48203E+03 

2.67165E+00 

3.89978E— 03 

— 9.04231E-07 

6.54370E— 11 

4.63902E— 15 

-9.63201E+03 

1.92568E+01 

6 

5.22642E+02 

3.50182E+00 

6.39254E— 03 

— 3.40202E-06 

8.93124E— 10 

— 9.30231E-14 

-4.24934E+03 

2.58749E+01 

7 

0 

4.96800E+00 

0 

0 

0 

0 

5.06219E+04 

— 9.14000E— 01 

8 

8.12266E+01 

4.87159E+00 

7.49130E— 05 

— 3.65807E — 08 

7.95819E— 12 

— 1.04187E— 16 

5.76456E+04 

1.09880E+01 

9 

1.94043E+03 

5.88669E-01 

7.14722E— 03 

— 2.90201E-06 

6.10537E— 10 

— 5.27109E— 14 

— 2.41961E+03 

4.44795E+01 

10 

— 3.78024E+02 

6.62079E-f 00 

3.29444E— 03 

-1.87396E-06 

5.15206E— 10 

— 5.54624E— 14 

2.17433E+04 

1.07487E+01 
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Table B2 lists coefficients B n j that appear in 
equation (35) for the six equilibrium constants K p j 
that correspond to reactions I through VI (eqs. (25) 


through (30)). The numerical values of B n j listed 
correspond to units for K p j that are nondimensional 
for j = 1, 3, and 6 and are in atmospheres for j = 2, 
4, and 5. 

Table B2. Equilibrium Constant Coefficients 


— 

B n ,j for " = 


j 

1 

2 

3 

4 

5 

6 

1 

5.09440E+03 

— 5.98336E+00 

1.65031E— 03 

— 4.45924E— 07 

4.89037E— 11 

8.42816E — 16 

2 

6.80400E+04 

-2.02038E+01 

— 2.04607E — 03 

1.83617E— 06 

— 6.39945E— 10 

8.33955E— 14 

3 

— 5.22510E+04 

1.04345E+01 

3.07030E-03 

— 1.40641E— 06 

3.75299E-10 

— 4.25248E— 14 

4 

-5.97702E+04 

1.26601E+01 

3.26442E— 03 

— 1.76622E— 06 

5.10028E— 10 

— 5.95591E— 14 

5 

-9.40571E+03 

3.89944E+00 

3.06797E— 05 

— 2.40321E— 07 

1.01593E— 10 

— 1.35585E— 14 

6 

— 2.10521E+04 

8.66294E— 01 

2.64558E— 03 

— 1.57780E— 06 

4.66330E— 10 

— 5.46540E — 14 


Table B3 lists coefficients D nk that appear in 
equation (63) for the viscosity of the seven species 
with k = 1, 2, . . . , 6, 10. The units for based on 
this set of coefficients are newtons-second per meter 2 . 
Table B4 lists coefficients E n ^ k in equation (66) for 
the thermal conductivity of the seven species with 
k = 1, 2, 6, 10. The units for X k based on 


this set of coefficients are joules per second-meter- 
kelvin. Note that the viscosity and thermal conduc- 
tivity data are only required in the region where wa- 
ter droplet growth occurs, which is generally at a 
temperature less than 500 K. Also, at these temper- 
atures the species corresponding to k = 7, 8, and 9 
do not appear in significant amounts. 


Table B3. Coefficients in Expression for Viscosity D nk 




D n,k 

for n = 


k 

1 

2 

3 

4 

1 

— 1.5371E— 05 

1.5429E— 07 

— 3.1981E— 10 

2.7958E — 13 

2 

— 2.5212E-07 

5.5024E-08 

— 1.2381E— 11 

— 5.5279E— 15 

3 

— 4.9877E-07 

7.8689E-08 

-6.9794E-11 

3.6438E — 14 

4 

— 5.7709E— 07 

8.8852E— 08 

— 7.0382E— 11 

3.4305E— 14 

5 

1.6225E-06 

2.9338E— 08 

— 1.9878E— 11 

1.1391E- 14 

6 

— 6.0917E— 08 

7.6311E — 08 

— 6.4962E— 11 

3.2374E — 14 

10 

— 8.4886E-07 

8.5390E— 08 

— 7.1595E— 11 

3.5702E — 14 


Table B4. Coefficients in Expression for Thermal Conductivity E n>k 


k 

B n ,k for " = 

1 

2 

3 

4 

1 

2 

3 

4 

5 

6 
10 

— 2.0450E— 01 

— 1.4403E— 05 
5.3712E— 04 

— 3.8505E-03 
— 5.5274E— 02 
— 3.0706E— 03 

— 1.8772E— 03 

1.7566E-03 
2.5166E— 05 
8.8425E— 05 
1.3168E-04 
1.2443E-03 
1.3215E-04 
1.1087E— 04 

— 4.6279E— 06 
1.3405E— 07 

— 1.5431E— 08 

— 1.2385E— 07 
-1.9310E-06 

— 1.4900E— 07 
— 6.9710E— 08 

4.2102E— 09 

— 1.1069E— 10 

— 1.8246E — 11 
8.1162E— 11 
1.3762E— 09 
1.0324E— 10 
2.7900E-11 
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Appendix C 

Computer Program 

This appendix presents a computer program, 
FIRACON, which has been developed to solve the 
finite-rate condensation problem with the solution 
technique of the present paper. The program is 


written in FORTRAN Version 5 language for the 
CDC® CYBER 180 series digital computer system, 
model 860, with Network Operating System. The 
program requires approximately 134000 octal loca- 
tions of core language and a typical case requires 
approximately 200 central processing unit seconds. 
Program input requirements and operating instruc- 
tions are summarized in the program. 


1 

2 
3 
A 

5 

6 

7 

8 
9 

10 

n 

12 

13 

1A 

15 
lb 
17 

16 

19 

20 
21 
22 
23 
2A 

25 

26 

27 

28 

29 

30 

31 

32 

33 
3 A 

35 

36 

37 
36 
39 
AO 


CONDENSATION 


PROGRAM F I RAC0N (INPUT, OUTPUT, TAPE 5-1 NPUT,TAPE6-0UTPUT, 

1 SAVDAT,TAPE7»SAVDAT,RDI5T,TAPE33«RDIST) 

C********* ****************** ******** ****** ***************************** 

C 

C FINITE RATE 

C 

C*******************************************************^^^^^^^^ 

c 

PROGRAM FIRACON PROVIDES A NUMERICAL SOLUTION FOR THE 
ADIABATIC EXPANSION OF COMBUSTION PRQ0UCTS WITH NUCLEATIUN 
AND WATER DROPLET GROWTH. THE NUMERICAL SOLUTION IS 
DIVIDED INTO ThREE successive parts, the first part 
INCLUDES THE COMPUTATION OF THE STAGNATION CONDITIONS 
IN THE COMBUSTION CHAMBER AND THE DETERMINATION UF THE 
MASS FLDW RATE, THE SECOND PART INVOLVES THE COMPUTATION 
OF ISENTROPIC FLOW FROM THE NOZZLE THROAT TO A POINT IN 
THE NOZZLE WHERE NUCLEATION JUST BEGINS. THE THIRD PART 
DEALS WITH THE COMPUTATION OF THE NUCLEATION PROCESS 
TO DETERMINE THE NUMBER UF DROPLETS FORMED AND THE 
SUBSEQUENT GROWTH OF THE LIQUID WATER DROPLETS, 


COMMON BLOCKS ARE USED EXTENSIVELY TO PASS INFORMATION 
BETWEEN THE VARIOUS SUBPROGRAMS. THE CONTENTS OF THE 
PRIMARY COMMON BLOCKS (PARTICULARLY AS THEY RELATE TO 
INPUT DATA AND OPERATING INSTRUCTIONS) ARE DESCRIBED 
IN THE DATA INPUT AND I N I T I AL I Z A T IUN SUBROUTINE, DATA 2 1 . 


THE USE 
TABLE: 


OF CUMHON BLOCKS IS SUMMARIZED IN THE FOLLOWING 


SUBPROGRAM 
FIRACON 
A V S X 
CCP 
CDEN 
CECOMP 
CEDATa 
CHG 


ABCDEFGHI 


COMMON BLOCK 
J K L M N 0 P 


QR5TUVWXY 


X 

X 

X X 


20 




41 

C 


CS X 

X 




X 



42 

C 


DATAil X X 

X XX 

X 

X 



X 


43 

C 


DBLIT X 

X 


X 


X 


X 

44 

C 


ELEMBAL 

X 


X 

X 

X 



45 

C 


EXPP 








4b 

C 


FITLAM 








47 

c 


FITMU 








48 

c 


F TEMP X 

X 

X 


X 



X 

49 

c 


GRL X 

X X 






X 

50 

c 


GRLPRP X X 

X X 

X 

X X 


X 

X 

X 

51 

c 


GRS X 

X X 




X 


X 

52 

c 


KPJ 








53 

c 


M I X LAM 





X 



54 

c 


Mixmj 

X 




X 



55 

c 


NUCRAT X 

X 




X 



56 

c 


PARTI X X 

X 

X X 

X 


X 

X 

X 

57 

c 


PART2 X X X 

X 

X X 

XX X 


X 

X X 

X 

56 

c 


PART3 X X XX 

X X 

X X 

X X 


X 

X 

X 

59 

c 


PLTOUT 

X 

X 

X XX 



X 

X 

60 

c 


PRTOUT X XXX 


X X 

X X 


X 

X 

X 

61 

c 


SATTEM X 





X 



62 

0 


XVSA X X 

X 







63 

c 










64 

c 


WHERE THE COLUMN HEADINGS REFER 

TO THE FOLLOWING 



65 

c 


common blocks* 








66 

c 










67 

c 


A) /ADJUST/ H) 

/C0NST5/ 

N J 

/K1K2ETC/ 


T) 

/ S IOCOM / 

6b 

c 


B) /ARECOM/ I) 

/EQCGM/ 

0) 

/OUT 1/ 


U) 

/SIGMAS/ 

69 

c 


C> / CMOuT 1/ J) 

/EXPDAT / 

P) 

/0UT2/ 


V) 

/svxz/ 

70 

c 


D) /CM0UT2 / K) 

/FEED/ 

OJ 

/PLTBLK/ 


W) 

/SWITCH/ 

71 

c 


EJ /CMOuT 3/ L) 

/GAMOUT/ 

R) 

/PLTOUT/ 


XJ 

/ TANDP / 

72 

c 


F) /CQEFCU/ M ) 

/ I NP T / 

S) 

/RESLTS/ 


YJ 


// 

73 

c 


G) /CONST/ 








74 

c 










75 



C0MM0N/5hITCH/SWEND,SWP0, SwGO 







76 



INTEGER S WEND * SWPO* S WGO 







77 



COMMON DUMMY ( 100U2 ) 








76 

c 










79 



CALL PSEUDO 








80 



CALL LEROY 








81 



CALL CALPLT(4.5*1.5*-3» 







82 


10 

CONT INUE 








P 3 



CALL D AT A 1 1 








64 

* 


IF ( 5 WEND « EQ • 1 ) GO TO 20 







85 



CALL PARTI 








86 



CALL PART2 








87 



CALL P ART 3 








88 



IF{ SWGO. EQ.l ) CALL PLTOUT 







89 



GO TO 10 








90 


20 

CONTINUE 








91 



CALL CALPLT<0. >0.*999) 








92 



STOP 








93 



END 
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SUBROUTINE AVSX(X,A> 

C 

C SUBROUTINE AVSX CALCULATES THE NOZZLE AREA AS A 

C FUNCTION OF X POSITION, 

C ACOEF CONTAINS THE SPLINE COEFFICIENTS DETERMINED 

C BY THE LIBRARY ROUTINE CSD5 

C 

COMMON/ ARECOM/N A, XAt ACOEF 

COMMON/ CONS T/TC,K,CAPR, Ml ,W (10) *RHQL»RBAR,PI, AL PH AC 
REAL K,M1 

0 I MEN S ION ACOEF ( 10,4 ) 

DIMENSION X A t 10 ) 

C 

NPD - NA-1 

XXI - X 

DO 10 I J* 1 , N P D 

IF C XA { I Jl ,LE.XXI.AND.XXI.LE.XA(IJ*in GO TO 20 
10 CONTINUE 
IJ - NPD 

20 H - XXI - XA(IJ) 

RBARX - ( ( ACDEF ( I J » 4) +H ♦ ACOEF ( 1 J * 3 ) )*H 
1 ♦ ACOEFUJ, 2))+H ♦ ACOfcFCIJ»lJ 
A - PI*RBARX**2 
RETURN 
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END 
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SUBROUTINE CCP l T , C P * C P 1 , CP C J 
C 

C SUBROUTINE CCP CALCULATES THE SPECIFIC HEAT CAPACITY FOR WATER 

C VAPOR AND THE CARRIER CAS 

C 

C0MMON/CON5T/TC ,K ,CAPR,M1 ,W( 10) , RHOL , R B A R , P I , ALPHAC 
REAL K * M 1 

CDMMON/COEFCO/COEF(10,8 i2) 

C0MM0N/5IGC0M /SIGMA (10),SIG,PSUP0 
C 

VAL - 0. 

J * 1 

IF ( T.LT. 1000. ) GO TO 10 
J * 2 

10 CONTINUE 
C 

C EQUATION (36) 

C 

DO 20 1*1,10 

VAL * VAL ♦ SIGMA (I)*(CO£F ( I » 1 » J ) /T ♦ C0EF(I,2,J) 

1 ♦ COE F ( 1 , 3 , J )+T ♦ CQEF( 1,4, j)*T**2 ♦ COEF (I , 5 , J ) *T**3 

2 ♦ COEFCI *6* JMT**4) 

IFCI.EO.il CPI - VAL 

20 CONTINUE 

CPC - VAL - CPI 
C 

C EQUATION 138) 

C 

CPI - tCPl/SIGHA(lM*4.184/W( 1) 

CPC • (CPC/CSIG - SIGMA(l) ))*4. 184 
C 

C EQUATION (24) 

C 

WC - 0. 

DO 30 1*2,10 

WC - WC ♦ S IGMAU )*W( I ) 

30 CONTINUE 

WC • WC/CS1G - 5IGMAUW 
C 

C EQUATION (40) 

C 

CPC - CPC/WC 
CP - 4 • 184* V A L 
RETURN 
END 


22 



SUBROUTINE CDE N UHO * S * T > 

SUBROUTINE CDEN CALCULATES THE MASS OENSITY OF 
THE GASEOUS MIXTURE BY INVERTING EQUATION <AA) 

COMMON/ CO NST/TC*K»CAPR*Ml»W(10)»RHOLtRBAR»P I * ALP HA C 

REAL K,M1 

’ COMMON/COEF CO/COE F< IQ ,8 ,2 ) 


10 

C 



11 



VAL « 0. 

12 



RP - CAPR/A.18A 

13 



J - 1 

1A 



IF( T.LT.1000.) GO TO 10 

15 



J - 2 

16 


10 

CONTINUE 

17 

C 



16 

C 


EQUATION (A6> 

19 

C 



20 



DO 20 I-l»10 

21 



VAL ■ VAL ♦ S IGMAl I>* (-CQEF ( I * ! 

22 



1 ♦ COE F ( i 1 3* J ) *T ♦ COEFU»A,J)' 

23 



2 ♦ COEF ( I *6* J J *T**A/A. ♦ CQEF1 

2 A 


20 

CONTINUE 

25 



VAL * A.18A*VAL 

2b 

C 



27 

C 


EQUATION ( A5 ) 

26 

C 



29 



VAL - S - VAL 

30 



VAL * VAL ♦ CAPR*SIG*ALOG<CAPR 

31 



DO 30 I*l»10 

32 



IFtSIGMAt n.EQ.O.) GO TO 30 

33 



VAL * VAL ♦ CAPR*SIGMA(I)*ALOG 

3 A 


30 

CONTINUE 

35 



VAL - -VAL/( CAPRASIGI 

36 



RHO - EXPP(VAL) 

37 



RETURN 

36 



END 


RP)*ALUG(T) 
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SUBROUTINE CECOMP 


1 



z 

C 




3 

C 

SUBROUTINE CECOMP CALCULATES THE CHEMICAL EQUILIBRIUM 



4 

C 

COMPOSITION OF COMBUSTION PRODUCTS 



5 

C 




6 


DIMENSION KP( 6 J » BB ( 6, 6 ) 



7 


REAL KP,K1,K2,K3,K4,K5,K6 



8 


REAL KPJ 



9 


COMMON/ K1K2ETC/K1,K 2 ,K3,K4,K5 ,K 6 



10 


COMMON/SIGMAS/SHiSOtSNfSC 



11 


COMM ON /CONS TS/F ACTOR i EPS 



1 2 


COMMON/RESLTS/S( 11) 



13 


COMMON/ INPT/TO,PO,XLA5T,PHU,RHC,RNO 



1 * 


CDMMON/TANDP/TXX.PXX 



15 


COMMON/SIGCOM/SIGMA (10) ,SIG,PSUPO 



16 

C 




17 

C 

EQUILIBRIUM CONSTANT COEFFICIENTS, B(N,J) 



Id 

C 

(N*l f 2f • 6 ) »( J-l»2 f . •• 96 ) ... TABLE 82 



19 

c 


\ 


20 


DATA BB/5.0943976E+3*-5.9833561,1.65Q3132E-3» 



21 


1 -4 . 4 592 4 15E-7 * 4. B903658E-11 , 0 . 4 281 563 E-16 »6 • 0O4OO26E ♦ 4, 



22 


2 -2. 020382 7E+1 »-2 . 04 6 065 2E -3 , 1. 8361667E-6, - 6 . 399 4 53 2E-10, 


- 

23 


3 8.33954B3E-14,-5.2250958E*4, 1 . 0 4 34 5 37E* 1 , 3 . Q7Q3042E-3, 


1 

24 


4 -1.4064069E-6,3.7529859E-10,-4.2524815E-14, 



25 


5 -5. 97 7015 2E* 4, 1.266 0142 E+ It 3. 264420 bE-3,-1 .7662154E-6 f 


j 

26 


6 5 . 100 2 83 5 E- 10 ,-5 .955 9096E-14 ,-9. 40 5 7111EO, 



27 


7 3.8994369, 3. 0679650E-5, -2. 40 3206 9E-7, 1.Q159289E-10, 



26 


8 -1. 3558463 E-14, -2 . 105 2 085 E* 4 , 8 . 6629408 E-l « 2 . 6 45 5B2 IE-3 , 



29 


9 -1 .5777998E-fe,4.6632962E-10,-5.4654014E-14/ 


] 

30 


DATA P5UP0/1 .01325E+5/ 

* 

; 

31 


DATA FACTOR/ 0.5/ 



32 


DATA EPS/1. E- 8 / 



33 

c 




34 

c 



- 

35 

c 

EQUATION (35) 



36 

c 




37 


DO 10 J * 1 , 6 


a 

38 


KP( J) * KPJ (BB( 1 , J) ,TXX ) 


- 

39 


10 CONTINUE 

: 


40 


K1 - l./KPU) 

- 


4 1 


K2 - (PXX/P5UPO)*KP(2) 



42 


K2 - 1./K2 



43 


K4 - KPC 3) / (PXX/PSUPO) 



4 


K 5 - KP(4)/(PXX/PSUP0) 

1 

= 

45 


K 3 - KP { 5 ) 



46 


K6 - KP ( 6 ) 


1 

47 


RCO - PH 1 1 / ( 2 . 0 ♦ 0 • 5*RHC ) 



4 fa 


SC « 1000.0/lRHC ♦ 12.0 ♦ ( 14 • 0*RN0 + 16.0J/RC0) 



49 


5H - RHC *S C 



50 


SO - SC/RCO 



- 51 


5 N « RNO*S 0 



52 


CALL ELEMBAL 

_ 

= 

53 


DO 20 1-1,10 



54 


SIGMAt I) - S(I) 

=■ 


55 


20 CONTINUE 



56 


5 IG • S(ll) 

= 

” 

57 


RETURN 

- 


58 


END 
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BLOCK DATA CEOATA 
C 

C HEAT CAPACITY COEFFICIENTS A(N,K) (N- 1 , 2 1 . . . , 6 ) 

C AND CONSTANTS OF INTEGRATION FOR ENTHALPY A(7,K> 

C AND ENTROPY A(8,K) ... TABLE B1 

C THESE COEFFICIENTS CORRESPOND TO THE SPECIES 

C H20, C 02 1 CO, 02 1 H2 , NZ, H, 0, OH , AND NO 

C IN THAT ORDER. 

C 

COMMON/ CO EFCO/COEF (10*8.21 

DATA < (COEFU.J.l), 1*1,10), J*1,B) 

1 /64 . 100011 *104 .75008,-1 .604 9999E+2 * 4 9, 4000 13 * 

2 -426.9,-151.54998,0. ,106.20001,165.30001,-9.1051247* 

3 7. 5915145, 3.3757382,9. 0628437, 6. 9432404,9. 8353258* 

4 8.8407937,4.968,5.0144162,6.1838936 ,8.4921514, 

5 -8 .4162903E-4, 2. 3892444E-Z ,-9 .9107277E-3 ,-3 .596 0249 E-3 , 

6 -7 .22 94 1 IE- 3, -8. 2 936 122 E-3, 0. ,-8 . 54 2 827 3E-4 , 2 . 60 374 78 E-3 , 

7 -9, 7741398E-3 , 6 .6 87 49 8 4 E-6 , -2 • 68 7 77E-5 , 1. 97942 88E-5 , 1 • 5 16025 5E-5 , 

8 8.4 80010 7E— 6,1.52779 5 6E— 5,0.,1 .74 32 IE -6,-5 .4702O64E-6, 

9 2 .320870 5E-5, -5. 195802 7E-9, 1.710 0074E-8 ,-l .495 7458 E-B , 

1 -1 .5427738E-8, 

2 -4. 452 796 3E -9 ,-l. 0301 117E-8 ,0. *-l .4656 169E-9 , 5. 6828683 E-9, 

♦ — 1.9844914E-08, 

3 1. 5576918E— 12,-4. 695 5142 E — 12 *4, 10 2563 7E— 12 ,5 . 2 051 27 7E-12 , 

4 1. 0 128202 E- 12, 2. 52 2 4 35 4E— 12 ,0 . , 4 .5512793E- 13 ,-2 .02564 14E-12 , 

5 6. 05 05 3 71 E-12 ,-60435.8 5, -96512 .03,-27911 .82,-2297.92,-244.86, 

6 -1519.28,50621 .9 4,5 7478 .92 , 645 4. 33 , 19 365 . 80 , 

7 2.0641,26. 1229,-2 .7t> 3 2, 10. 134 7 ,-24.446, 

0 -3.2292, -C. 9140, 10.4430,0. 62 13, 3*9782/ 

DATA (( COEF( I, J,2) , 1-1,10) , J-1,8) 

1 /3491. 9931 ,-2226 . 6 402 ,-3 . 50 34 1 44 E+2 ,— 1 5 39 .2227, 

2 148 2.02 60,522 .64 184 ,0 . , 8 1 . 2 26626 , 1940 , 4 33 9 ,-3 . 7802350E* 2, 

3 -4. 64 1B407, 13. 739656, 6. 3 17705 8, 10. 738 034, 2. 67 1645 3, 

4 3.5018215, 4. 96 8, 4. 87 15 058 ,0.50866853,6.6207915, 

5 1.688102 2E-2,2*5 53989 2 E-3 , 3 . 2 644 3E-0 3 ♦ -1 . 3 869361 E-3 , 

6 3.8997751E-3,b.3925364E-3,0. , 7 .49 1304 3E-5 , 7 . 1472198E-3 , 

+ 3.2944365 E-3, 

7 -7 . 39 5 3 04 8E -6,-1. 425 0656E-6, -1. 6662 24 7E-6, 5.811 5599E-7, 

8 -9.0 423126-7,-3.402 02 27E-6 , 0 . ,-3 . 65 806 7E-8 *-2 . 902 OOoE-b , 

+ - 1. 87 3961 IE-6 , 

9 1.6850683E-9,3.7628479E-10,4.1344159E-10,-5.58070l8E-llf 

4 6.54370046—11,8 .9312 3 96E-10 ,0 . , 7. 9581929E- 12 ,6 .1053683E-10 , 

1 5. 1520591E-10, 

2 -1.5 8 3412 E-13,- 3. 877904 9E-14, -4. 03865896-14,-3. 2845032E-15, 

3 4. 63 90234E-15, -9. 30230976-14,0 • *-l . 04 10673E-16 ,-5 .271O910E-14, 

4 -5. 5462447E -14, -77426. 11, -85 3 15. 28, -26246. 25, 5834. 96, -9632. 01, 

5 -4249.34,50621.94,57645.58,-2419.61,21743.30, 

6 77.4501 ,-34 .7509,9 .4698,-16.4089,19.2568, 

7 25.8749,-0.9140, 10.9860,44.4795,10.7487/ 

END 


ORIGINAL' PAGE IS 
OF POOR QUALITY 
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SUBROUTINE CHG1T ,HG ,HG1 ) 

C 

C SUBROUTINE CHG CALCULATES THE SPECIFIC ENTHALPY 

C OF THE GASEOUS MIXTURE AND OF WATER VAPOR 

C 

COMMON/COEFCO/COEF ( 10,8*2) 

COMMON/CONST/TC»K,CAPR,Ml,WtlO) ,RHOL,RBAR,P I,ALPHAC 
REAL K,M1 

COMMON/ S IGCOM/ SIGMA (10) , SIG,P5UP0 
C 

VAL * 0, 

J - 1 

1 F ( T * LT • 1000 • ) GO TO 10 
J - 2 

10 CONTINUE 


16 

17 

18 

19 

20 
21 
22 

23 

24 
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26 

27 

28 

29 

30 

31 


C 

C EQUATION i 42 ) 

C 

DO 20 1*1,10 

VAL * VAL ♦ SI GMA ( X)*{COEF(I,lfJ)+ALOG(T) ♦ COEF ( 1 ,2 » J ) *T 

1 * C0EFU,3,J)*T*A2/2. ♦ COE F ( 1 , 4, 4 ) *T**3/3 . 

2 ♦ C0EF(I,5,J)*T**4/4« ♦ CQEF (1 ,6, 4\+T**5/5. ♦ C0EFU«7»JII 
IF(I.EO,i) HG1 - VAL/SIGMAd) 

20 CONTINUE 

HG - 4.184*VAL 
C 

C EQUATION (43) 

C 

HG 1 - HG1*4 • 1B4/W ( 1 ) 

RETURN 

END 
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SUBROUTINE CS(S,RHQ,T) 

C 

C SUBROUTINE CS CALCULATES THE SPECIFIC ENTROPY 

C OF THE GASEOUS MIXTURE 

C 

COHMON/CONST/TCtK,CAPR,Ml,W(10) *RHOL,RBAft,PI, ALPHAC 
REAL K , M 1 

C OMMON/COEFCO/ COEF ( 10,8,2) 

COMMON/ S IGCOM/ S IGMA (10), SIG,P5UP0 
C 

VAL * 0, 

RP - CA PR/ 4 * 18 4 
J - 1 

IF ( T . L T. 1000* ) GO TO 10 
J - 2 

10 CONTINUE 
C 

C EQUATION (44) 

C 

DO 20 1-1,10 

VAL - VAL ♦ SIGMAC I >♦ (-COEF ( I , 1, J )/T ♦ (C0EF(I,2,J) - RP)AALQG(T) 

1 ♦ C0EF(I,3,J)*T ♦ C0EF(I,4, J)*T**2/2. ♦ COEF d * 5, J ) ATAA3/3. 

2 ♦ C0EF(I,6,J)+T*A4/4» ♦ C0EF(I,8,J)) 

20 CONTINUE 

VAL - 4 • 184* VAL 

VAL * VAL - CAP RA SI G* ALOG ( RHO) - C AP R* SI G* ALQG ( CAPR/PSU PQ ) 

DO 30 1-1,10 

1F(SIGMA( I) • EQ • 0* ) GO TO 30 

VAL - VAL - CAPRA SI GMA( I)*ALOG(SIGMA(I) ) 

30 CONTINUE 
S • VAL 
RETURN 
END 
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SUBROUTINE DATAII 


1 

2 

3 

4 

5 

6 
7 
e 

9 

10 

11 

12 

13 

14 

15 

16 
17 
16 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 
3 l 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 


C 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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SUBROUTINE DATAII HANDLES THE DATA INPUT AND 
INITIALIZATION, THE SUBROUTINE CSOS IS A PART OF THE 
MATHEMATICAL LIBRARY AT LANGLEY, FOR MORE INFORMATION 
SEE "MATHEMATICAL AND STATISTICAL SOFTWARE AT LANGLEY'S 
CENTRAL SCIENTIFIC COMPUTING COMPLEX DOCUMENT N2-3C. 


DIMENSION ACOEFI 10*4)*DY(10I ,XA ( 10) ,RWA <10 ) ,DELSTR<10) 
DIMENSION Y ( 10 I 
DIMENSION WK ( 79 ) 


"EXPDAT" VARIABLES 


NEXP - NUMBER OF EXPERIMENTAL DATA POINTS TO BE 
DISPLAYED ON PRESSURE PLOT, IF NEXP IS 
ZERO* NO EXPERIMENTAL DATA IS PLOTTED, 

IF NEXP IS NOT ZERO, THE EXPERIMENTAL 
DATA IS INCLUDEO ON PLOT TYPE 1. 

XEXP - X POSITIONS FOR WHICH EXPERIMENTAL DATA 
IS TO BE DISPLAYED ON PRESSURE PLOT 

PE XP - EXPERIMENTAL VALUES OF PRESSURE TO BE 
DISPLAYED ON PRESSURE PLOT 

PREF - REFERENCE PRESSURE - I.E., THE VALUE PEXP/PREF 
IS DRAWN ON THE PRESSURE PLOT 

C0MM0N/EXPDAT/NEXP,XEXP(50) , P EX P { 50 ) ,P RE F 

"SWITCH" VARIABLES 


S WEND - END OF DATA SWITCH 

0 - DATA ENCOUNTERED, PROCESS IT 

1 - END OF DATA ENCOUNTERED, STOP 

SWPO - PRINTED OUTPUT SWITCH 

0 - DO NOT PRINT OUTPUT 

1 - PRINT OUTPUT (DEFAULT) 

SWGO - GRAPHIC OUTPUT SWITCH 

0-00 NOT PLOT OUTPUT 
1 - PLOT OUTPUT (DEFAULT) 

COMMON/ S WITCH/ 5 WEND , SWPO, SWGO 
INTEGER 5WEND, SWPO, SWGO 

"ARECOM" VARIABLES 


NA - NUMBER OF X POSITIONS DEFINING NOZZLE GEOMETRY 

XA - NOZZLE X COORDINATES 

ACOEF - SPLINE COEFFICIENTS DEFINING AREA VS X 
COMMON/ ARECOM/NA,XA,ACOEF 
"FEED" VARIABLES 

HFUEL - ENTHALPY OF THE FUEL AT THE FEED 
TEMPERATURE, CAL/MOL 

TF EEO - TEMPERATURE OF THE FEED, K 

NATOM - NUMBER OF ATOMS OF CARBON IN A MOLECULE 
OF FUEL 

COMMON/FEED/HFUEL, TFEED *NATOM 
"PLTBLK " VARIABLES 


PLTBLK CONTAINS THE MINIMUM AND MAXIMUM VALUES 
FOR ALL PLOT OUTPUT, IF ANY PARAMETERS ARE NOT 
SPECIFIED, THE CORRESPONDING PLOT IS SCALED 
AUTOMATICALLY 

XMIN 
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DISTANCE ALONG THE NOZZLE 


83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 
117 
116 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


XMAX 
PHIN 

PMAX - NORMALIZED PRESSURE 
TM IN 

TM AX - NORMALIZED TEMPERATURE 

YM IN 

YMAX - MASS FRACTION OF LIQUID WATER 
DTMIN 

OTMAX - TEMPERATURE DIFFERENCE 
DNMIN 

ONMAX - NUCLEATION RATE 
RM IN 

RMAX - DROPLET RADIUS 
PRMIN 

P RMAX - PRESSURE RATIO 
TLMIN 

TLMAX - TEMPERATURE OF LIQUID 


RHOMIN 

RHOMAX - DENSITY 
SMIN 

SMAX - ENTROPY 
UMIN 

UMAX - FLOW VELOCITY 
XMMIN 

XMMAX - MACH NUMBER 

IPSLCT - PLOT SELECTION INDICATOR. IF THE VALUE 15 
ZERO FOR A GIVEN INDEX, THE PLOT IS NOT 
GENERATED. IF THE VALUE IS ONE FOR A GIVEN 
INDEX* THE PLOT IS GENERATED. THE INDICES 
FOR ALLOWED PLOTS ARE AS FOLLOWS* 

1 - NORMALIZED PRE55URE 

2 - PRESSURE RATIO 

3 - NORMALIZED TEMPERATURE 

4 - MASS FRACTION 

5 - TEMPERATURE DIFFERENCE 

6 - NUCLEATION RATE 

7 - CRITICAL RADIUS 

8 - DROPLET RADIUS (AND R+) 

9 - TEMPERATURE OF LIQUID 

10 - DENSITY 

11 - ENTROPY PRODUCTION 

12 - FLOW VELOCITY 

13 - MACH NUMBER 

14 - COMBINATION OF INDICES 1* 3* AND 12 

15 - COMBINATION OF INDICES 4* 5, AND 6 

COMM ON/P LTBLK/XM IN, XMAX, PM IN ,P MAX »TM IN » TM AX * 

1 YMIN*YMAX,DTMIN*OTMAX,DNMIN,DNMAX,RMINf RMAX* 

2 PRM IN ,PRMAX, TLMIN, TLMAX, RHOMIN* RHOMAX, 

3 SM IN, SMAX, UMIN, UMAX, XMMIN, XMMAX, IP SLCT (15 J 

"CONST" VARIABLES 


TC - CRITICAL TEMPERATURE OF WATER, K 

K - BOLTZMANN CONSTANT, J/MOLECULE-K 
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167 

C 




168 

C 

CAPR 

- 

UNIVERSAL GAS CONSTANT* J/NQL-K 

169 

c 




170 

c 

Ml 

- 

MOLECULAR MASS OF WATER* KG/HOLECULE 

171 

c 




172 

c 

w 

- 

MOLECULAR WEIGHT TABLE. KG/MOL 

17 3 

c 




174 

c 

RHOL 

- 

DENSITY OF WATER* KG/H**3 

175 

c 




176 

c 

R6AR 

- 

SPECIFIC GAS CONSTANT FOR WATER* R/Wl* J/KG-K 

177 

c 




178 

c 

PI 

- 

AREA OF UNIT CIRCLE 

179 

c 




180 

c 

ALPHAC 

- 

THERMAL ACCOMMODATION CUEFF I Cl ENT FOR CARRIER 

181 

c 



GAS INTERACTION WITH WATER DROPLETS 

182 

c 




183 


COMMON/ CO N5T/TC*K*CAPR*M1*W(10) * RHQU t RBAR, P I , ALPHAC 

184 


REAL K t Ml 



185 

c 




186 

c 

"ADJUST" 

VARIABLES 

187 

c 




188 

c 




189 

c 

GAMMA1 

- 

RATIO OF SPECIFIC HEATS 

190 

c 




191 

c 

QC 

- 

CONDENSATION COEFFICIENT 

192 

c 




193 

c 

ALPHA 

- 

CONSTANT IN EQUATION (78J 

194 

c 




195 

c 

BETA 

- 

LANGMUIR PARAMETER 

196 

c 




197 

c 

DSTAR 

- 

DIAMETER OF NOZZLE THROAT. M**3 

198 

c 




199 

c 

ASTAR 

- 

AREA OF NOZZLE THROAT* M**3 

200 

c 




201 

c 

DT2 

- 

TEMPERATURE STEP FOR PART 2* K 

202 

c 




203 

c 

JMIN 

- 

BEGINNING NUCLEATION RATE FOR PART 3* 

204 

c 



DROPLETS F0RMED/M*43-S 

205 

c 




206 

c 

DELX 

- 

STEP WIDTH FOR PART 3 CALCULATIONS* M 

207 

c 




208 

c 

EPS 

- 

CONVERGENCE ARRAY 

209 

c 

1 

- 

NEWTON ITERATION TO SOLVE EQUATION (571 

210 

c 

2 

- 

NEWTON ITERATION TO SOLVE EQUATION (111) 

211 

c 

3 

- 

PRESSURE CONVERGENCE IN PART 1 

212 

c 

4 

- 

MASS FLUX CONVERGENCE IN PART 1 

213 

c 

5 

- 

PRESSURE CONVERGENCE IN PART 2 

214 

c 

6 

- 

NEWTON ITERATION FOR T(J-l) IN PART 3 

215 

c 

7 

- 

TEMPERATURE CONVERGENCE* EQUATION (133)* IN PART 

216 

c 

6 

- 

NEWTON ITERATION FOR X AS A FUNCTION OF A 

217 

c 




216 

c 

I OUT 

- 

OUTPUT LEVEL 

219 

c 

0 

- 

STANDARD SUMMARY OUTPUT 

220 

c 

1 

- 

EXTENDED OUTPUT 

221 

c 

2 

- 

DEBUG OUTPUT 

222 

c 




223 

c 

DTI 

- 

TEMPERATURE STEP FOR PART 1* K 

224 

c 




225 

c 

CVRAT 

- 

CONVERGENCE RATIO FOR DELS/CV* EQUATION (126) 

226 

c 




227 

c 

ARAT 

- 

CONVERGENCE RATIO FOR A'/A, EQUATION (131) 

226 

c 




229 

c 

JSTOP 

- 

MAXIMUM NUMBER OF BANDS IN PART 3 

230 

c 



DIMENSION STATEMENTS LIMIT THIS VALUE TO 400. ANY 

231 

c 



SMALLER NUMBER WILL PROVIDE A MEANS TO STOP THE 

2 32 

c 



CALCULATIONS SHORT OF THE END OF THE NOZZLE. 

233 

c 




234 

c 

JDB 

• 

DEBUG PRINT CONTROL FOR PART 3 

235 

c 



IF THIS VALUE IS DEFINED* 'EXTENSIVE' PRINTOUT 

236 

c 



IS PROVIDED FOR DEBUGGING BOTH BEFORE AND AFTER 

237 

c 



THE DBLIT ROUTINE FOR ALL BAND NUMBERS GREATER 

236 

c 



THAN OR EQUAL TO JDB 

239 

c 




240 

c 

CP 

- 

SPECIFIC HEAT CAPACITY OF MIXTURE* J/KC-K 

241 

c 




242 

c 

CPI 

- 

SPECIFIC HEAT CAPACITY UF WATER* J/KG-K 

243 

c 




244 


C OMMON/AD J US T/GAMMA 1,0C» ALPHA f BET At DSTAR* AS TAR. DT2» 

245 


1 JMIN *DE LX i 

, EPS ( 8) .1 OUT* DTI. CVRAT* ARAT* JSTOP. JDB *CP.CP1 

246 


REAL JMIN 



247 

c 




246 

c 

" I NPT" 

VARIABLES 

249 

c 




250 

c 
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251 

C 

TO - INITIAL TEMPERATURE, K 

252 

C 


253 

C 

PO - INITIAL PRESSURE, N/M**2 

254 

c 


255 

c 

XL AST - FINAL X, M (USED TO TERMINATE CALCULATIONS) 

256 

c 


257 

c 

PHI 1 - EQUIVALENCE RATIO 

256 

c 


259 

c 

RHC - ELEMENTAL RATIO OF HYDROGEN TO CARBON 

260 

c 

IN MIXTURE 

261 

c 


262 

c 

RNO - ELEMENTAL RATIO OF NITROGEN TO OXYGEN 

263 

c 

IN MIXTURE 

26* 

c 


2 65 


COMMON/ INPT/TO,PO»XLAST »PHI1,RHC,RN0 

2 66 

c 


267 

c 

NAMELIST PARAMETERS 

268 

c 


269 

c 


270 

c 

DESCRIPTIONS OF MOST NAMELIST PARAMETERS CAN BE 

271 

c 

FOUND ELSEWHERE IN THIS SUBROUTINE. THE FOLLOWING 

2 72 

c 

PARAMETERS ARE LOCAL TO THIS ROUTINE. 

273 

c 


27* 

c 

IPLTRD - INTEGER FLAG TO INDICATE IF AUTO-SCALING 

275 

c 

IS TO BE USED ON THE PLOTS. 

276 

c 

0 - USE AUTO-SCALING 

277 

c 

1 - DO NOT USE AUTO-SCALING. THE NAMELIST 

278 

c 

"PLTOVR" MUST BE INCLUDED IN THE INPUT DATA. 

279 

c 


280 

c 

RNA - WALL COORDINATES FOR DEFINING TUNNEL GEOMETRY 

281 

c 


282 

c 

DELSTR - BOUNDARY LAYER DISPLACEMENT THICKNESSES FOR 

283 

c 

DEFINING TUNNEL GEOMETRY 

28* 

c 


265 


NAMELIST/DATAIN/TO,PQ,X LA ST, SWPOfSWGOiQC, ALPHA, BETA, 

286 


1 DSTAR,CVRATiARAT,0Tl,I0UT,DT2»JMIN,DELX,EPS, 

287 


2 JSTDP, JOB , PHI 1, RHC, RNO, IPLTRD ,NA,XA ,RW A, DELSTR,ALPHAC, 

28B 


3 NEXP,XEXP,PEXP,PR£f,IP5LCT»HFUEL,TFEED,NAT0M 

289 


NAMELIST/ PL TO VR/XM IN ,XMAX,PH IN, PM AX, TM 1N.TMAX, 

290 


1 YMIN,YMAX,OTMIN,DTMAX,DNMIN,DNMAX,RMIN,RMAX, 

291 


2 PRM IN,PRMAX,TLMI N, TLMAX , RHOMIN , RHOMAX , 

292 


3 5MIN,5MAX,UMIN,UMAX,XMM in,xmmax 

293 

c 


29* 

c 

SET CONSTANTS 

295 

c 


296 


DATA TC, K, CAPR, Ml *W,RH0L/6 *7.286, 1. 361E-2 3, 8. 31* ,2 .99E-26, 

297 


1 0.01 8, 0.0 **,0.028, 0.03 2, 0.00 2, 0.028, 0.001 ,0.016, 0.0 17, 

298 


1 0.030, 1.E3/ 

299 


OATA DY/10* 1. / 

300 

c 


301 

c 

SET DEFAULTS FOR ADJUSTABLE CONSTANTS 

302 

c 


303 


DATA ALPHAC/1./ 

30* 


DATA TO/O. / 

305 


DATA HFUEL,TFEED,NAT0M/-17895.,298.15,1/ 

306 


OATA I P5LCT/15*0/ 

307 


OATA NEXP/O/ 

308 


OATA XEXP/92., 116. ,1*0., 16 *.,188 .,212* ,22* . ,236. ,2*8 . , 

309 


1 260., 272. ,28*., 308 .,320 .,332. ,3** ., 356 . ,39 2 *26. , 

310 


2 *6*. ,500. ,536., 616. ,27*0./ 

311 


OATA OC, ALPHA, BETA ,DSTAR/ 1.0, 9.0, 

312 


1 2.0,0.01128*/ 

313 


DATA S WEND/O/ 

31* 


OATA I PLTRD/O/ 

315 


OATA XMIN,XNAX ,PMIN,PHAX»TMIN»THAX ,YMIN»YHAX, 

316 


1 DTMIN,DTMAX,DNMIN,0NNAX,RNIN,RMAX/l**O./ 

317 


DATA P RN I N, PRM AX, TLMIN, TLMAX, RHOMIN, RHOMAX, 

318 


1 UMIN,UMAX,SMIN,SMAX,XMMIN,XMMAX/12*0./ 

319 


DATA CVRAT/l.E-3/ 

320 


DATA ARAT/l.E-3/ 

321 


DATA EPS/8*l.E-6/ 

322 


DATA PH 11/0.8/ 

323 


OATA RHC/*./ 

32* 


DATA RNO/ 3.7619/ 

325 


DATA I OUT / 1/ 

326 


DATA DTI/1.0/ 

327 


DATA DT2/0.2/ 

326 


DATA NPAS5/0/ 

329 


DATA 5 WPQ , SMGO/ 1,1/ 

330 


DATA JST0P/100/ 

331 


DATA JOB / 1000 / 

332 

c 


333 


1 FORMAT(1H1,20X » 2 2HC0ND ENS AT ION EVOLUTION//) 

33* 

c 
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335 

336 
J37 

338 

339 
3*0 
3*1 
3*2 
3*3 
3** 
3*5- 
3*6 
3*7 
3*8 
3*9 

350 

351 

352 

353 
35* 

355 

356 

357 

358 

359 

360 

361 

362 

363 
36* 


SET COMPUTED CONSTANTS 

RBAR ■ CAPR/Wll) 

PI - AC0$<-1.) 

READ (5tDATAlNtEND«10) 


10 

IF { EOF (5 ) 

• NE, 

,0) 

GO 

TO 


WRITEI33) 

TO, 

>P0 




1F( IPLTRD 

• NE, 

.0) 

READ (5 

20 

I F ( E OF ( 5 ) 

• EO ■ 

.0) 

GO 

TO 

30 

SWEND - 1 





*0 

CONTINUE 






I F ( TO , NE • 

0 , ) 

GO 

TO 

50 


TO * 1500* 

CALL F TEMP ( TFL AME ) 

NTENS * ITFLAHE - 95.J/1Q. 

TO - 10.*FL0AT(NTENS> 

50 CONTINUE 

IF (NPASS.EQ.O) RRI TE ( 6 1 1 ) 

IF INPASS ,EQ, 0) WRI TE (6 ♦DATAIN) 

NPAS5 - NPASS ♦ 1 

I ERR - 0 

NPM - 10 

S - 0,0001 

IPT - -1 

DO 60 I *1 * N A 

YU) - RWA(I) - DELSTRII) 

60 CONTINUE 

CALL C5DS(NPM*NA*XA*Y»DYiS* IPT * ACOEF • WK « IE R R ) 

RETURN 

END 
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SUBROUTINE DBLIT 


1 

2 

3 

4 

5 

6 
7 
B 
9 

10 

11 

12 

13 

14 

15 

16 
17 
lb 

19 

20 
21 
22 

23 

24 

25 

26 
27 
26 

29 

30 

31 

32 

33 

34 

35 

36 

37 
36 

39 

40 


C 

C 

C 

C 

C 

C 


c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 


SUBROUTINE DBLIT DETERMINES THE TEMPERATURE FOR WHICH 
TWO INDEPENDENT CALCULATIONS OF ENTROPY ARE EQUAL AND 
THE DENSITY FOR WHICH THE MASS FLOW AND CROSS 
SECTIONAL AREA ARE CONSISTENT 

COMMON/S IGCOM /SIGMA (10),$ I G,P5UP0 

COMMON I, J ,A( 400) ,X(400) ,DELN(400) , CAP JX ( 400 ) * ft (400) » 

1 MU( 400) ,TL(400) ,T(400> ,P ( 400 ) , RST AR ( 4 00 ) , 

2 U( 400) ,S(400) ,SL(4 00) ,HL(400) ,RHO(4QO) ,CAP JY1400J * 

3 F(400) *G(400)»YS( 400 ) ,DELY(400) >DELS(400) » 

4 RPERM(400) *TS(400) ,0UTMCH(400) ,CAPPJ(400) 

REAL MU 

COMMON/CQN$T/TC,K,CAPR ,M 1 «W (10 ) , RHOL * RBAR * P I » ALPHAC 
REAL K »M 1 

COMMON/ AO JUST/GAMMAl.QC, ALPHA, BETA, DSTAR,A STAR *DT2* 

1 JMIN*DELX,EPS(8)» IOUTtDTl*CVRAT.ARAT» jSTOPtJDBtCPtCPi 
REAL JMIN 

COMMON/ DU T 1/50 »H0, MOOT *RHO STRUTS TAR 
REAL MDOT 


S GPP - (S(J + 1) - SL(J*1))/(1. - YS ( J ♦ 1 ) ) 
10 CONTINUE 

EQUATION (44) 


CALL CS(SGP*RH0U*1) *TCJ*11) 

CALL CCP(T(J*1) ,CP*CP1,DUH) 

EQUATION (126) 

U50CV - (SGPP - S GP ) / ( CP - CAPR*SIG) 
IFiABS (D50CV) .LT.CVRAT) GO TO 20 

ADJUST T, EQUATION (127)* AND GO BACK 

T C J ♦ 1 1 - T< Jtl )+EXPP(D50CV) 

GO TO 10 


41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 
67 
6b 

69 

70 

71 

72 

73 

74 


C TEMPERATURE ITERATION HAS CONVERGED 

C 

20 CONTINUE 

CALL CHG ( T ( J ♦ 1 ) *HGP*OUM) 

C 

C CALCULATE SPECIFIC ENTHALPY USING EQUATION (128) 

C 

HP - (1. - YS(J*1))*HGP ♦ HLUU) 

C 

C CALCULATE FLOW VELOCITY USING EQUATION (129) 

C 

U( J*l) • SQRT(2.*(H0 - HP) i 
C 

C CALCULATE CORRESPONDING AREA USING EQUATION (130) 

C 

AP * MDOT* ( 1 * - YS(J + l))/( RHQ (J + l)*U(J + l) ) 

APRAT - AB 5 ( AP - A(J + l))/AU + l) 

IF(APRAT.LT.ARAT) GO TO 30 
C 

C ADJUST RHO* EQUATION (132), AND GO BACK 

C 

RHQ ( J ♦ 1 ) - MOOT* ( 1 • - YS ( J ♦! ) )/(A(J+l)*U(J+l)) 

GO TO 10 
C 

C DENSITY ITERATION HAS CONVERGED 

C 

30 CONTINUE 
C 

C CALCULATE PRESSURE USING EQUATION (47) 

C 

P(J+1> • RHO( J*1)*CAPR*SIG*T(J*1) 

CALL $ATTEM(TS( J+1),P( J+l) ) 

RETURN 

END 



1 

z 

3 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 


SUBROUTINE ELEMBAL 
REAL K1»K2»K3,K4»K5,K6 
C0HM0N/RESLTS/SU1) 

COMMON/SIGMA5/SH,SO,5N,SC 
COMMON /CONSTS /FACTOR, EPS 
COMMON/K1K2ETC/K1,K2,K3,K4,K5,K6 
1 FORMAT ( 1 *♦* ELEMBAL EXCEEDS 50 ITERATIONS 
S ( 2 ) -SC *F ACTOR 
S ( 1 ) -SH/2 • 0 

S (4>-( SO-SC-SI 1>-S(2) J/2.0 
S (6J-SN/2 ,0 

5(11 }-(SH*S0*SN45C-5( 1 J-5(2) )/2.0 
NITER-0 
10 CONTINUE 
Y CLD- S ( 2 ) 

ZOLD - S ( 4 ) 

NITER-NITER+1 
IF(NITER.LT.51) GO TO 20 
WRITE (6t 1) 

STOP 

C 

C COMPUTE THE REST OF THE VARIABLES ....... . 

C 

20 CONTINUE 

S(3)-S(2)*SQRTIK2*S(11)/S(4)J 
5(5)-5(l)*S(3J/ ( S ( 2 ) *K 1 } 

5(10) - SQRT(K6*S( 4 )*5 (6 ) ) 

5(6) - (SN-S( 10))/2.0 
5(9)-SQRT(K34S(4 )*$(5) ) 

$(7)-SQRT(K4*S(5)*S(ll) ) 

S(8)*SQRT(K5*S(4)*S(11M 
S (1) -l SH-Z.0*5t 51-5 (7 1-5(9) 1/2.0 

sm-sc-sm 

5(4)»(S(4) + (50“SC-5(l)-5(2)-5(8)-S(9)-S(10) )/2.0) /2.0 
S(li)-SU)*S(2>*Sj3l*S(4)*5C5)*S(6MSC7)*5(8)»5(9MS(10> 
C 

C CHECK FOR CONVERGENCE 

C 

DY-Y OLD-5 ( 2 ) 

OZ - ZOLD - S ( 4 ) 

I F ( NITER .EQ. 1) GO TO 10 
IF ( AB5 ( DY/S ( 2 J ) .LE. EPS) GO TO 30 
GO TO 10 
30 CONTINUE 

IFUBS IDZ )/S(4) .LE.EPS ) GO TO 40 
GO TO 10 
40 CONTINUE 

5(3)-S(2)45QRT(K2*S( ll)/5(4) ) 

S(5)-5(1)*S(3)/(S(2)*K1) 

S( 10) - 5QRT(K6*S(4)*S(6) ) 

S(6) - ( SN— 5 ( 10 ) ) /2 . 0 
5(9)»5QRT(K3*S(4)*S(5) ) 

S< 7)-5QRT( K4*S (5>*S( 11 ) J 
5(6)-SQRT(K5*5(4)*S(U)) 

RETURN 

END 


ORIGINAL FAGS IS 
OF POOR QUALITY 



FUNCTION EXPPIVAU 


1 

2 

3 

4 

5 

6 
7 

e 

9 

10 

n 


c 

C FUNCTION EXPP IS A UTILITY USED TO AVOID THE UNDERFLOW 

C MESSAGE ASSOCIATED WITH VERY SHALL ARGUHENTS 

C 

EXPP - 0* 

IF(VAL.LT.-670.) GO TO 10 
EXPP - EXP(VAL) 

10 CONTINUE 
RETURN 
END 


" 

1 


SUBROUTINE F ITLAM ( L AKBD A , T ) 



2 

C 




3 

C 

SUBROUTINE FITLAH DETERMINES THE THERMAL 

CONDUCT 

* 


C 

OF EACH SPECIES AS A POLYNOMIAL IN TERMS 

OF TEMPI 

- 

5 

C 



“ 

6 


dimension lambdaiioj 



7 


REAL LAMBDA 


- 

8 


DIMENSION C OEF { 4 1 10 ) 



9 

C 




10 

C 

COEFFICIENTS IN EXPRESSION FOR THERMAL 


= 

11 

c 

CONDUCTIVITY EtN,K) ... TABLE 84 



12 

c 



’ 

13 


DATA COEF/- 0.2045,1. 7566E-3*-4 .6279E-6,4 .2102E-9, 


14 


1 -1.4403E-5,2.5166t-5*1.3405E-7*-1.1069E-10* 



15 


2 5.3712E-4*B.8425E-5*-1.543lE-8,-l .8246E-11, 



16 


3 -3.8 505 E-3 *1 .3 168E-4, -1.2 38 5E -7,8 .1162E-11 * 


* 

17 


4 -5. 52 74E-2,1.2443 E-3 »-1.9 3 ICE-6,1 .3762E-9, 


- 

18 


5 -3.0706E-3* 1. 321 5E -4 * - 1 . 490UE-7 . 1.0 32 4E-10, 


* 

19 


6 12 + 0. t - 1 .8 7 72 E-3 * 1 • 1087E-4 ,-6 . 971 QE-8 , 2 .790GE-11/ 


20 

c 




21 

c 

EQUATION (661 



22 

c 




23 


DO 20 1-1*10 



24 


LAMBDA ( I ) - 0. 



25 


DO 10 J- 1 * 4 



26 


LAMBDA(I) - LAMBDA(I) ♦ C OEF ( J , I ) + T* + ( J-i ) 



27 


10 CONTINUE 



28 


20 CONTINUE 



29 


RETURN 


_ 

30 


END 



— 

I 


SUBROUTINE F1TMU(MU,T> 


2 

C 



3 

C 

SUBROUTINE F1TMU DETERMINES THE VISCOSITY OF EACH 


4 

C 

SPECIES AS A POLYNOMIAL IN TERMS OF TEMPERATURE 


5 

C 



6 


DIMENSION MU ( 1 0 ) 


7 


REAL MU 


8 


DIMENSION COEF { 4 » 10 ) 


9 

C 



10 

C 

COEFFICIENTS IN EXPRESSION FOR VISCOSITY 


11 

c 

D(N*K) ... TABLE B3 


12 

C 



1 3 


DATA C0EF/-1. 5371E-5 *1 .5429E-7 *-3. 19 81E -10 , 2 . 7958E-13 * 


14 


1 -2.5212E-7,5.502 4E-8*-1.2381E-U*-5*52 79E-15* 


15 


2 -<». 98? 7E-7, 7, 8689 E-8 *-6.9 79 4E- 11,3.64386-1 4, 

I 

16 


3 -5. 7709E-7, 8 . 8852E-8 ,-7.0 382E-11, 3.43Q5E-14, 

~£ 

17 


4 1.6225E-6 *2 .9338E-8 *-1.96786-11* 1. 1391 E-14 , 


16 


5 -6.0917E-8,7.6311E-b*-6.4962E-ll*3.2374E-14, 


19 


6 12+0.,-8.4886E-7*8.5390E-8,-7.1595E-ll*3.5702E-14/ 


20 

c 



21 

C 

EQUATION (63) 


22 

C 


' 

23 


DO 20 1-1*10 


24 


MU ( I ) - 0. 


25 


DO 10 J - 1* 4 


26 


MU(I) - MU(I) ♦ COEF ( J * I ) + T++ ( J- 1 J 


27 


10 CONTINUE 


26 


20 CONTINUE 

-- 

29 


RETURN 


30 


END 


34 
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57 
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C 

C 

C 

C 

C 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


SUBROUTINE F TEMP ( TFL AME ) 

SUBROUTINE FTEMP COMPUTES THE FLAME TEMPERATURE OF 
HYDROCARBON COMBUSTION IN AIR WITH OR WITHOUT 
OXYGEN ENRICHMENT ••••••■ 

CONSTANT PRESSURE F 0 RMUL AT ION . • . . 

REFERENCE TEMPERATURE IS 296,15 KELVIN 


ORIGINAL PAGE IS 

of; poor quality 


COMMON/ INPT/TO, PO , XLA 5 T , PHI 1 , RHC,RNO 
COMMON/TANDP/TXX, PXX 
COMMON/RESLTS/ S ( 1 1 ) 

COMMON /FEED/HF UEL, TFEED iNATOM 
COMM ON /COEF CO/ COEF (10,0 ,2) 

1 FORMAT!/, 2X, 'DID NOT CONVERGE IN 20 ITERATIONS') 

2 FORMAT! /,2X, 'ENTHALPY OF THE FEED (CAL/KG.) - SF12.4) 

3 FORMAT!/, 2X, 'FlAME TEMPERATURE (KELVIN) « ',F8.2) 

4 FORMAT* /,2X, 'DENSITY OF THE MI XTURE ( KG/CU.M )•* ,E12 ,6 ) 

5 F0RMAT(/t2X, 'ENTHALPY OF THE MIXTURE (CAL/KG.)- *, F12.4) 
VAL(I,J,T) - (COEF ( I »1,J ) * ALOG ( T) + COEF ( I ,2, J ) *T 

1 ♦ COEF ( I ,3, J)*T4*2/2. ♦ COEF( I,4,J)*T**3/3. 

2 ♦ COEF (I ,5, J) *T**4/4. ♦ COEF ( I , b , J ) *T * *5 /5 • ♦ C0EF(I,7,J)) 

HFUEL - ENTHALPY OF THE FUEL AT THE FEED TEMPERATURE (CAL/MOL) 
TFEED - TEMPERATURE OF THE FEED (KELVIN) 

N - NUMBER OF ATUMS OF CARBON IN A MOLECULE UF FUEL, 

TFO * FIRST GUESS OF FLAME TEMPERATURE* KELVIN. 

PO - PRESSURE OF THE FEED (BARS) 

COMPUTE SIGMA(H)* SIGMA(O), SIGHA(N)* AND SIGMA(C) 

IN 1 KG OF THE MIXTURE 


EQUATION (8) 

SO ■ 500 , 0/ ( 8,0+7,0*RNO+PHI 1*1 12 , 0 + RHC ) / t 4,0+RHC) ) 
EOUATIGN (9) 

SC - 2.04PHI 1*S0/ ( 4.0 + RHC) • 

EQUATION (11) 

SH-RHC *S C 

EQOATION (10) 

5N-S0*RN0 

ENTHALPY OF THE FEED.... 

FUEL HAS THE FORMULA C*N H*M 

THERE ARE N MOLS OF CARBON IN ONE MOL OF THE FUEL, 
HFEEO-HFUEL * SC * FLOAT ( N ATOM ) 

HFEED-HFEED+VAL (4*1* TFEED )*50/2.0*VAL (6,1, T FEE D ) ♦ SN/ 2 . 0 

COMPUTE THE EQUILIBRIUM COMPOSITION AT TFO AND PO.. 
COMPUTE THE ENTHALPY 

T X X - TO 
PXX - PO 
CALL CECOMP 
TFO - TXX 

CALL CHG(TFO, HO, DUMMY) 

NOTE: SUBR CHG RETURNS H IN JOULES/KG. 

T 1 -TFO- 100 • 0 

DH1-HFEED-H0/4.154 

DO 10 ITER-1,20 

TXX - T1 

CALL CECOMP 

CALL CHG( Tl, HI, DUMMY) 

DH2-HFEED-H1/4. 184 

TFLAME-(DH2*TF0-DH1*T1 )/ (DH2-DH1) 

IFIABS (TFLAME-T1) .LE. 0,001) GO TO 20 
TF0-T1 
Tl-TFLAME 
DH1-0H2 
10 CONTINUE 
WRI TE ( 6, 1 ) 

STOP 



B3 

20 CONTINUE 


84 

OENCTY-PO/8 

• 314/St 1D/TFLAME 

85 

wRITE(6,2) 

HFEED 

8b 

WR 1 TE ( 6» 3 ) 

TFLAME 

07 

*RITE(6,4) 

DENCTY 

88 

WRITE (6*5) 

HI 

89 

RETURN 


90 

END 



“ 

1 


SUBROUTINE GRL(ZZ) 


2 

C 


“ 

3 

c 

SUBROUTINE GRL CALCULATES THE GROWTH RATE 


4 

c 

FOR LARGE DROPLETS (I.E., Z > 1.1) 

— 

5 

c 



6 


COMMON/ EQCOM/DELTAT*CA PL |NU*P$I tDROT » 


7 


1 TR* LT ILDA, KNt THETA *LAMBDA,PRiX I> OMEGA tCAPL AM 


6 


REAL NU f LTILDA,KN, LAMBDA 

- 

9 


COMMON/CON$T/TCtK,CAPR»Hl ,w(10) tRHOL » RBAR tP I t ALPHAC 

= 

10 


REAL K » H 1 

- 

ll 


COMM ON /AD J UST/ 0 AMMA 1 »Q C * ALPHA » BETA fDSTARt AS TAR fOTZ* 

“ 

12 


l JMIN»DELX,EPS t8) t IOUT,OTliCVRAT,ARAT» JSTOP tJDBfCPiCP 

- 

13 


REAL JM IN 


14 


COMMON I» JtAl400) ,X ( 400 ) , DELNt 400 ) ,C AP JX ( 400) »R<400) « 


15 


1 MU(400) »TLC400) ,T(400) ,P(400) »RSTARt400) » 

“= 

16 


2 Ut 400) f S (400) t SL (400) t HL(400) t RHO ( 400 ) i C AP JY ( 400) » 


17 


3 F < 4 00 ) # G(400> tYS( 400) t DEL Y( 400) * DELS (400) , 


16 


4 RP ERM 1 400 ) *T$(4Q0> tOUTMCH ( 400 ) , CAPP J ( 400 ) 


19 


REAL MU 


20 

c 



21 


ZOLO • ZZ 


22 


RAT - R( I ) /RSTARt J) 


23 


CAPDEL - CAPLAM4DELTAT 


24 


ZPRIME - (THETA**3/(THETA ♦ 1. > ) ♦ ALOGt RAT ♦ THETA) 


25 


1 ♦ (1. /(THETA ♦ 1.) ♦ OMEGA )*ALOG ( RAT - 1.) 


26 


2 ♦ (OMEGA ♦ 1. - THETA ) *RAT ♦ 0.5*RAT+*2 ♦ CAPOEL 

— 

27 


10 CONTINUE 


26 

c 



29 

c 

EOUATION (112) 


30 

c 



31 


FZ2 - (THETA**3/(THETA ♦ 1 • ) ) *ALOG ( ZOLD ♦ THETA) 


32 


1 ♦ (1. /(THETA ♦ 1.) ♦ 0MEGA)4AL0G(Z0LD - 1.) 

■iz 

33 


2 ♦ (OMEGA ♦ 1. - THETA ) *ZOLD ♦ 0.5*Z0LD**2 - ZPRIME 


34 

c 



35 

c 

EQUATION (113) 

— 

36 

c 



37 


FZ2P « THETA**3/ ( (THETA ♦ l.)*(ZOLD ♦ THETA)) 


38 


1 ♦ (1. ♦ OMEGA* ( THETA ♦ l.))/((THETA ♦ l.)*(ZOLD - 1. 


39 


2 ♦ ZOLD ♦ OMEGA ♦ 1. - THETA 

- 

40 

c 


= 

41 

c 

EOUATION (114) 


42 

c 



43 


ZNEW * ZOLD - FZ2/FZ2P 

— 1 

44 


I F ( A BS ( ZNEM - ZOLD) /ZNE W. LT .EPS ( 2 ) ) GO TO 20 

- 

45 


ZOLD - ZNEW 

= 

46 


GO TO 10 



47 


20 CONTINUE 

— 

46 


ZZ - ZNEW 


49 


RETURN 


50 


END 



36 


■ill 
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SUBROUTINE GRLPRP 
C 

C SUBROUTINE GRLPRP CALCULATES SEVERAL VARIABLES 

C PREPARATORY FOR EVALUATION OF THE GROWTH RATE FOR 

C LARGE DROPLETS 

C 

COMMON/S IGCOH /SIGMA ( 10) ,$IG,PSUPO 
DIMENSION WORK ( 10 J 

COMMON/EQCOM/DELTAT,CAPL,NU,PSI,DRDT , 

1 Tft* LT ILD A , KN» THET A , LAMBD A » P R , XI , OME G A , CAP LAM 
REAL NU,LTILDA,KN, LAMBDA 
COMMON/I NPT/TO,PO,XLAST, PHI 1,RHC,RN0 
COMMON/ S M I TCH/SW£ND,SWP0,5WGG 
INTEGER S WEND » SWPQ » 5 wGQ 

COMMON/ 0UT2 /U2 , RH02 * CA P J2 ,T2 ,P 2 ,RS T AR2 ,X2 v A2, S IG2U1) ,OUNC 
C0MM0N/CM0UT3/NENT 

COMMON/ CUT 1/ SO, HQ, MOOT, RHQSTR,T STAR 
REAL MOOT 
REAL KNST 

COMMON/ AD JUST/GAMMA1,QC, ALPHA, BET A, DSTAR* AS TAR, DT2, 

1 JMIN,DELX,EPS(BJ , IOUT ,DT 1 , C VRAT, ARATt J5T0P t JDB ,CP , CP l 
REAL JMIN 

C0MM0N/C0NST/TC,K,CAPR,M1 ,W ( 10 J ,RHOL , RB AR , P I , ALPHAC 
REAL K , Ml 

COMMON I , J , A (400) «X (400) ,DELN (400) ,CAP JX (400) ,R(400) » 

1 MU ( 400) *TL ( 400 ) ,T(400) , P ( 400 ) , RSTAR ( 400 ) , 

2 U1400) ,5(400) ,SL(400) ,HL(400) ,RH0(4 00) ,CAP JY(400) , 

3 F ( 400 ) ,G(400) ,YS(400) , DEL Y (400) , DELS (400) , 

4 RPERM(400J ,TS< 400) , OUTM CH ( 400 ) , C APP J(400) 

REAL MU 

REAL MUU 
C 

CALL 5ATTEM(TSP,PU)) 

CALL CHG (TIJ ) ,HGX,HG1) 

C 

C EOUATION (56) 

C 

CAPL - HG1 - 4,2E+3*T( J) ♦ 17 . 1 1 75 3658 E *6 
CALL CCP(T( J) , CP, CPI, CPC) 

WAVG - 0. 

WC - 0. 

00 10 IJ-1,10 

WAVG - WAVG ♦ W( IJ) *SIGMA( IJ) 

IF(IJ.NE.l) WC - WC ♦ W(I J)*$IGMA(IJ) 

10 CONTINUE 

WAVG - WAVG/SIG 

WC - WC/(5IG - SIGMA(l) ) 

WRAT * W ( 1 ) /WAVG 
Y 1 - S I GMA ( 1 ) / S IG 
C 

C EQUATION (39) 

C 

GAMMA1 - CPi/(CPl - CAPR/Wd)) 

C 

C EQUATION (41) 

C 

CAMMAC - CPC/(CPC - CAPR/WC) 

C 

C EQUATION (74) 

C 

FF - Y 1* SORT ( WRAT ) ♦ (1* - Y1 ) *SuRT ( WC /WAVG ) * ( GAMMA1/GAMMAC ) 

1 4 ( ( GAMM AC ♦ 1 • )/ ( GAMMA1 ♦ 1. )) * (CPC/CPI) *ALPHAC 
C 

C EQUATION (78) 

C 

NU - ((CAPR/W(1))*TSP/CAPL)4(ALPHA ♦ 0.5 - ((2. - 00/(2.400) 

1 4UGAMMA1 ♦ 1. )/(2.4GAMMAl) )4( (CP14TSP)/CAPL) 

2 4(FF/Y1)/SQRT(WRAT>) 

TR - T ( J )/TC 

CALL F I TMU ( WORK , T ( J ) ) 

CALL MIXMU(MUU,WORK) 

C 

C EQUATION (68) 

C 

LTILDA - 1.54MUU*SQRT(CAPR*S1G4T( J) )/P( J) 

C 

C EQUATION (71) 

C 

KNST « LTILDA/(2.*RSTAR( J)) 

C 

C EQUATION (75) 

C 


37 



83 


THETA - 2.*BETA*KN 

BA 


CALL FITLAMl WORK,T 

85 


CALL M IXLAM (LAMBDA 

86 

C 


87 

C 

EQUATION (69) 

88 

C 


89 


PR - CP1*MUU/LAM8D 

90 

C 


91 

c 

EQUATION (77) 

92 

c 


93 


XI - S QRT ( 8 • *P I ) * ( 

9 A 

c 


95 

c 

EQUATION (73) 

96 

c 


97 


OMEGA - XI*I 1. - N 

96 

c 


99 

c 

EQUATION (76) 

100 

c 


101 


CAPLAH - L AHBDA + 1 T 

102 


RETURN 

103 


END 


1 

z 

3 

A 

5 

6 

7 
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9 
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16 

17 
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23 

2 A 

25 

26 

27 

28 

29 

30 

31 

32 

33 
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35 
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37 

38 
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A2 
A3 
AA 
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C 

C 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


FUNCTION GRS(ZZ> 

FUNCTION GRS CALCULATES THE GROWTH RATE FOR 
SMALL DROPLETS ,U.E.» 2 < 1.1) 

COMMON/S IGCOH/ SIGMA (10),SlG,PSUP0 
C0 MMON/EQCOM/DELTaT,CAPL,NU,PSI»DROT, 

1 TR,LTILDA,KN,THETA,LAMBDA,PR,XI,0MEGA,CAPLAH 
REAL NU»LTILDA,KN, LAMBDA 

COMMON/ CON ST/TC »K , CAP R , M 1 , M ( 10 ) , RHQL , RBA R» P I , ALPHAC 
REAL K f Ml 

COMMON/ AO JUST /GAM MAI ,GC, ALPHA ,8 ETA* D STAR ,A5TAR,DTZ , 

1 JM IN, DE LX, EPS I 8) » I OUT, DTI *C VRAT,ARAT» JS TOP , JDB • CP »CP 1 
REAL JM IN _ _ 

COMMON I , J , A ( A 00 ) ,X ( AOO ) , DE LN ( A 00 ) , C AP JX ( A 00 J ,R(AOO), 

1 MU (AOO) ,TL( AOO) ,T( AOO) »P (AOO) ,RSTAR(AOO) , 

2 U( AOO) , S ( AOO ) ,SL(AOO) ,HL(AOO) , RHO ( AGO) ,CAPJY(AOO), 

3 F(AOO) , G ( AOO ) ,YS(AOO),DELY(mOO),DELS(AOO), 

A RPERM { AOO) ,TS(AOO) ,UUTMCH ( AOO ) ,CAPP J ( AOO) 

REAL MU 

CALL CCPtTl J) , CP, CPI, CPC) 
nAVG * 0. 

WC - 0. 

00 10 IJ«1,10 

WAVG * WAVG ♦ W( I J)*SIGMA< I J 1 
IF(IJ.NE.I) WC - WC ♦ W ( I J ) *S IGMA ( I J ) 

10 CONTINUE 

WAVG - WAVG/SIG 

WC - WC/( SIG - SlGMA(l) ) 

WRAT » W ( 1 ) / w A VG 
Y1 » SIGMA(li/SIG 

EQUATION (39) 

GAMMA1 - CP1/(CP1 - CAPR/wm) 

EQUATION (Al) 

GAMMAC - CPC/(CPC - CAPR/WC) 


EQUATION ( 7A ) 

FF » Y1*SQRT l W RAT ) ♦ (l. - Y1 ) *SURT ( WC/WAVG )* (GAMHA1/GAMMAC) 

1 + { (GAMMAC ♦ 1 • ) / ( GAMM Al «• 1 .))♦( CPC/CP 1 ) *ALPHAC 

EQUATION (80) 

PSI -FF*P(J)/SQRT(2.*P1*CAPR*5IG*T(J))*( (GAMMA 1 ♦ 1 • ) / ( 2# *GAMHA l ) ) 
1 *CP1*RSTAR(J)/(CAPL*RH0L*( 1. - NU ) ) * ( TS ( J ) - TUJ) 

PSIBAR - PSI/(R(I)*RSTAR( J)) 

ORDT - ( RSTAK ( J*l) - RSTAR( J) ) /D6LTAT 

EQUATION (79) 

SOL - U ♦ RSTAR( J ) /RS T AR U* 1 > * ( 2 Z - 1 . ) *E XPP (PS I B AR*DELTAT ) 

1 - ( l./( PSIBAR*RSTAR( J ♦ 1 ) ) ) *ORDT* ( E XPP ( PS I B AR*DELTAT) - 1»> 

GRS - SOL 

RETURN 

END 


38 
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REAL FUNCTION KPJ(B,T) 

DIMENSION B ( to ) 

C 

C 

C EVALUATE EQUATION (35) FOR THE EQUILIBRIUM CONSTANTS 

C 

KPJ - EXPP(B(1)/T ♦ B(2> ♦ BI3)*T ♦ 6(4)*T**2 
1 ♦ B ( 5 ) +T++ 3 ♦ B (6 ) *T+*4 ) 

RETURN 

END 
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SUBROUT INE M IX LAM ( L AM BOA , LAMBDA I ) 

C 

C SUBROUTINE MIXLAM CALCULATES THE THERMAL CONDUCTIVITY 

C OF THE MIXTURE 

C 

DIMENSION LAMBDAK6) 

REAL LAMBDA * LAMBD A I 

COMMON/S IGCON/SIGMAt 10) ,SIG,PSUPO 

C 

C 

C EQUATION (67) 

C 

SUM - 0. 

DO 10 1-1,10 

IFtLAMBDAI ( D.EQ.O. ) GO TO 10 
SUM - SUM ♦ SIGMA(I)/LAMBOAI(I) 

10 CONTINUE 

LAMBDA « 5IG/SUM 
SUM - 0. 

DO 20 1-1,10 

SUM - SUM ♦ S IGMA ( I ) +L AMBDAI ( I ) 

20 CONTINUE 

LAMBDA - 0.5MSUM/SIG ♦ LAMBDA) 

RETURN 

END 
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SUBROUTINE M IXMU ( MU, MU I ) 

C 

C SUBROUTINE MIXMU CALCULATES THE VISCOSITY OF THE 

C GASEOUS MIXTURE 

C 

COMMON / CONS T/ T C, K ,C APR, Ml, W ( 10) * RHOL • R BAR , P I , ALPHAC 
REAL K • Ml 
DIMENSION MU I ( 10 ) 

REAL MU ,MUI 

COMMON/S IGCOM/SIGMA < 10 ) ,$IG,P5UP0 
C 

c 

c EOUATION (64) 

C 

MU - 0, 

00 20 1-1,10 
SUM - 0. 

DO 10 J-1,10 

IF(J.EQ.I) GO TO 10 

I F ( MUI ( J ) • EQ • 0« ) GO TO 10 

PHI - (1. ♦ SORTCMUHI )/MUI( J) ) *(W ( J)/H(1M *+0«2S)**2 
PHI - P H I / ( SORT ( 8* ) ♦SQRKl* ♦ M(I)/W(J))) 

SUM - SUM ♦ SIGMA (J ) *PHl 
10 CONTINUE 

IF(SIGMA(I)«E0*0« ) GO TO 20 

MU - HU ♦ MU1U)/(1. ♦ SUM/ SIGMA (I)) 

20 CONTINUE 
RETURN 
END 
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16 
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26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 
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47 

48 

49 

50 

51 

52 


C 

C 

C 


C 

C 

C 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


SUBROUTINE NUCRAT (CAP J ,T ,P ) 

SUBROUTINE NUCRAT CALCULATES THE NUCLEATION RATE, J 
COMMON/S IGCOM/SIGMA(10),SIG,PSUPO 

COMMON/CONST/TC»K,CAPR ,M1 ,W (10) ,RHOL, RBAR , P I ,ALPHAC 
REAL K , Ml 

COMMON/ AD JUS T/GAM MAI, OC, ALPHA, BE TA,DSTAR, AS TAR, 0T2, 

1 JMIN,DELX,EPS18)» IQuT ,0T 1 ,C VRA T, AR AT , JSTOP ,JDB ,CP,CP 1 
REAL J MIN 


EQUATION (57) 

PINF - E XP P ( 55 • 897 - 6641. 7/T - 4*4864*AL0G(T)> 

PBAR - (P/PINF)*(SIGMA(1)/SIG) 

TR - T/TC 

EQUATION (58) 

SIGMAX - (82.27 ♦ 75.612*TR - 256. 889*TR**2 ♦ 95. 928+TR443 > 
1 M.E-3 

EQUATION (59) 

RSTARX - 2.*SIGMAX/{RH0L*RBAR*T*AL0G(PBAR) ) 

CALL CHG( T ,HGX , HG1 ) 

EQUATION (56) 

CA PL • HG1 - 4.2E*3*T ♦ 17 . 1 175 3o58E +6 
RHOX - P/ ( CAPRAS IG + T) 

GAMMA 1 - CP1/(CP1 - CAPR/Wll)) 

EOUATION (62) 

0 « 2 . * ( GA MMA 1 * 1. )*CAPLMCAPL/(RBAK*T) - 0.5)/(RBAR*T 
1 * ( GAMM A1 ♦ 1.)) 

WAVG - 0* 

DO 10 1*1,10 

WAVG * WAVG + S IGMA( I )*W ( I ) 

10 CONTINUE 

WAVG - WAVG/5IG 
mR AT - ( W ( 1 ) /WAVG ) 

EQUATION (61) 

CAPJ - (l./U. ♦ Q))*QC*WRAT4SQRT(2.4SIGMAX/(PI*M1**3)) 

1 *( RH0X**2/RH0L)*EXPP(-4.*PI*S IGMAX*R$TARX**2/(3.*K*T1 ) 

CAPJ • CAP J MS IGM A ( 1 ) /5 IG ) **2 

RETURN 

END 
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OF POOR QI 

1 



SUBROUTINE PARTI ^ 

2 

C 



3 

c 


STAGNATION CONDITIONS ANO NASS FLOW RATE 

4 

c 



5 



CON NON/CO NS T/TC»KiCAPRfMl » W ( 10 ) » RHOL , RB AR , P I , ALP HAC 

6 



REAL K » Ml 

7 



CONNON/GANOUT/GANBAR 

b 



COMMON/ AD JUST/GAMNA1,GC* ALPHA, BET A,D$TARi AS TAR f DT2* 

9 



L JhIN,DELX,EP$m. IOUT,DTliCVRAT,ARATt JSTOPi JOB » CP t CPI 

10 



REAL JMIN 

11 



COMMON /SWITCH/S WEND *SWPO,SwGO 

12 



INTEGER SWENDtSWPO*SWGO 

13 



COMMON/ INPT/TOy PO »XLAST • PHI 1 « RHCt RNQ 

14 



COMMDN/TANDP/TXX.PXX 

15 



C0MH0N/CUT1/S0 , HQ» MDQT ♦ RHOSTR » T 5TA R 

16 



REAL MOOT 

17 



COMMON/CMOUT1/NENTRY ,DT ,RHOO» T1 , HI ,U 1, US TAR »URH0*RH0F1 

16 



COMMON/ S IGCOM/ SIGMA ( 1 0 ) , S IG » P S UPO 

19 

c 



20 

c 


COMPUTE ASTAR, AREA AT NOZZLE THROAT 

21 

c 



22 



NENTRY - 0 

23 



ASTAR - PI/ 4.*QSTAR+*2 

24 

c 



25 

c 


COMPUTE HO « HG USING EQUATION (42) WITH T • TO 

26 

c 



27 



T XX - TO 

28 



PXX • PO 

29 



CALL CECOMP 

30 



CALL CHG( T0,H0,DUM1 J 

31 



RHOO * P0/(CAPR*SIG*T0) 

32 

c 



33 

c 


COMPUTE SO - SG USING EQUATION (44) WITH T - TO* RHO 

34 

c 



35 



CALL C5 ( 5 0 * RHOO * TO I 

36 

c 



37 

c 


MAKE INITIAL GUESS FOR T 

38 

c 



39 



T 1 - 0 • 9*T 0 

40 



T XX - T 1 

41 



PXX « PO 

42 



DT « DTI 

43 



URHOX - 0. 

44 


10 

CONTINUE 

45 



NENTRY * NENTRY ♦ 1 

46 


20 

CONTINUE 

47 



CALL CECOMP 

48 

c 



49 

c 


COMPUTE DENSITY USING EQUATION (45) 

50 

c 



51 



CALL CDEN(RHOFltSOtTl) 

52 

c 



53 

c 


COMPUTE PRESSURE USING EQUATION (47) 

54 

c 



55 



P - RH0F1*CAPR*SIG*T1 

56 



DIF « AB5 ( PXX - P) 

57 

c 



58 

c 


TEST FOR CONVERGENCE ON PRESSURE 

59 

c 



60 



IF(DIF/P.LT. EPSOM GO TO 30 

61 



PXX » P 

62 



GO TO 20 

63 


30 

CONTINUE 

64 



CALL CHGf TL «HGF«DUM1) 

65 

c 



66 

c 


COMPUTE FLOW VELOCITY USING EQUATION (103) 

67 

c 



68 



UF • SORT ( 2 • * ( HO - HGF ) ) 

69 



HI « HGF 

70 



U1 - UF 

71 

c 



72 

c 


COMPUTE MASS FLUX 

73 

c 



74 



URHO - RHOF 1 *UF 

75 



GAMBAR - AL0G(T1/T0)/AL0G(P/P0> 

76 



GAMBAR - l./(l. - GAMBAR) 

77 



IFdOUT.GT.UAND.SWPa.EO.il CALL PRTQUT(U2> 

78 

c 



79 

c 


IF PAST A MAXIMUM FOR MASS FLUX BRANCH TO CHANGE 

80 

c 


RESOLUTION TO LOCATE THE MAXIMUM 

81 

c 



82 



IF ( URHO. LT. URHOX ) GO TO 40 


RHOO 


41 



83 

84 

85 

86 

87 

88 

89 

90 

91 
9 Z 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 


URHOX - URHO 
II - Tl- DTI 
T XX • Tl 
GO TO 10 
C 

C PAST A MAXIMUM* RESTORE LAST T AND DECREASE 

C DELTA T - GO BACK AND TRY AGAIN 

C 

40 TST - ABSIURHO - URHQXI/URHO 
I F ( TST • LT * E PS ( 4 M GO TO 50 
Tl - Tl ♦ 2.0*DT1 
TXX - Tl 
0T1 - 0* 5 *QT1 
URHOX - 0* 

GO TO 10 
C 

C HAVE FOUND A MAXIMUM MASS FLUX 

C SET TSTAR - T* RHOSTR • RHD* USTAR - U 

C 

50 TSTAR - Tl 

RHOSTR - RH0F1 
USTAR - UF 
C 

C COMPUTE MASS FLOW RATE, MOOT USING EQUATION 1104) 

C 

MOOT * RHOF l*UF*AST AR 
IF(SWPO.EQ.l) CALL PRTDUTUtO) 

RETURN 

END 
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SUBROUTINE PART2 

ISENTRQPIC EXPANSION FROM NOZZLE THROAT 

COMMON/5 IGCOM/SICMA (10 ),$IG,P SUP 0 
COMMON/ 5VX2/X22»Y22 
COMMON/GAMOUT/GAMBAR 
COMHON/SWITCH/SWEND,SWPO»SWGU 

COMMON/P LTOUT/NQUT, TOUT (600) »POUT(600> ,X0UT(600> « 


10 



1 OELTI (600>,UOUT ( 600 ) ,0UTMACt600J 

11 



INTEGER SWEND,SWPO,SWGO 

12 



COMMON/ INPT/TO.PO, XL AST, PHI 1,RHC*RN0 

13 



COMMON/TANDP/TXX,PXX 

1 A 



C OMM ON /OUT 1/SO, HO, MOOT, RHOSTR,T STAR 

15 



REAL HOOT 

16 



COMMON/ CMOUTi/NENTRY,DT,RHQO,Tl,Hl,Ul,U STAR, URHOtRHOFl 

17 



COMMON/ AD JUST /GAMMA1, OC, ALPHA, B ETA, DSTAR, AS TAR, DT2, 

18 



1 JMIN,DELX,EPS(B),I OUT ,DT1,CVRAT,ARAT,JST0P,JDB,CP,CP1 

19 



REAL JMIN 

20 



C0MM0N/C0N$T/TC,K,CAPR,H1 ,H ( 10 ) , RHOL , RBAR, P I ,ALPHAC 

21 



REAL K »M 1 

22 



COMMON/ 0UT2 /U2 , RH02 , CAPJ2,T2,P2,RSTAA2,X2,A2,SIC2(ll) , OUMC 

23 



COMMON/CM 0UT2/NE,T2N,SIGMAX,U2W,RH0,P, PI NF,PBAR,TR,H2W,RSTAR, 

2 A 



1 Q,J 

25 



REAL J 

26 

C 



27 


1 

FORMAT* 1H0, '**♦** ARRAY SIZE EXCEEDED IN PART 2 **♦♦*') 

28 


2 

FORMAT (1H0, END OF NOZZLE ENCOUNTERED WITH INSUFFICIENT 

29 



1 ,* VALUE OF J ♦****• ) 

30 

C 



31 

c 


START WITH T - TSTAR AND RHO - RHOSTR 

32 

c 



33 



NENTRY - 0 

3 A 



NENTRY2 * 0 

35 



IBR » 0 

36 



XOLD - 0. 

37 



T2W » TSTAR 

38 

c 



39 

c 


DECREASE T BY DT2 

AO 

c 



A 1 



NE • 0 

A2 



P - RHOSTR* CAP R*S IG*TSTAR 

A3 


10 

CONTINUE 

AA 



NE - NE ♦ l 

A5 



T2N - T2W - DT2 

A6 



TXX - T2W 

A7 


20 

CONTINUE 

A8 



PXX - P 

A9 



CALL CECOMP 

50 

c 



51 

c 


CALCULATE DENSITY USING EQUATION (45) 

52 

c 



53 



CALL CDENCRH0ZZ,S0,T2W) 

5 A 

c 



55 

c 


CALCULATE PRE5URE USING EQUATION (47) 

56 

c 



57 



PZZ - RHQZZ*CAPR*5IG*T2W 

5b 



DIF - ABS ( P - PZZ) 

59 



IF(DIF/PZZ.LT.EPS(5)> CO TO 30 

60 



P • PZZ 

61 



GO TO 20 

62 


30 

CONTINUE 

63 

c 



6A 

c 


CONVERGENCE IN PRESSURE ITERATION 

65 

c 


CALCULATE ENTHALPY USING EQUATION (42) 

66 

c 



67 



CALL CHG (T2 W ,HG , OUM ) 

68 

c 



69 

c 


CALCULATE FLOW VELOCITY USING EQUATION (103) 

70 

c 



71 



U2W - SORT ( 2 * * ( HO - HG ) ) 

72 



H2W - HG 

73 



RHO - RHO ZZ 

7 A 



NOUT - NE 

75 



TOUT ( NOUT ) - T2W/T0 

76 



POUT (NOUT I - P/PO 

77 



UOUT ( NOUT ) - U2W 

78 



CALL SATTEM ( TS , P ) 

79 



DELTT ( NOUT ) - TS - T2W 

80 



IFtOELTT(NOUT) • IT . 0. ) DELTT(NOUT) - 0. 

B 1 

c 



82 

c 


CALCULATE CPI USING EQUATION (38) 


43 



83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 
117 

ue 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 
13b 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 
16 3 

164 

165 

166 


C 

CALL CCP(T2W,CPtCPl*CPC> 

C 

C CALCULATE GAHHA1 USING EQUATION (39) 

C 

GAMHA1 * CP1/(CP1 - C APR/W ( 1 ) ) 

OUTHAC(NOUT) - GAHMA1*C APR4S 1 G*T2W 
OUTH AC ( NOUT > - U2W/$QRT< OUTMAC ( NOUT )) 

C 

C CALCULATE PINF USING EQUATION (57) 

C 

PINF » EXP P ( 55 *897 - 6641.7/T2W - 4 .4 8644 ALOG ( T2W U 
PBAR - (P/PINF)*(SIGMA(1)/SIG) 

C 

C CALCULATE AREA USING EQUATION (105) 

C 

ALST - MD0T/(RH0*U2H) 

C 

C GET X POSITION CORRESPONDING TO THIS AREA 

C 

CALL XVSA(ALSTtXLST) 

IF ( ( XLS T - XOLO ) . GT. 0. 1 ) DTZ - 0.5*DTZ 
XOlD - XLST 

I F ( NOUT • GT* 595 ) NRITEUtl) 

IFiN0UT.GT.595) STOP *100 MANY POINTS* 

XOUT(NOUT) - XLST 
X2 - XLST 

GAHBAR » AL0G(T2N/T0)/AL0G(P/P0) 

GAMBAR - l./(l. - GAHBAR) 

NENTRY • NENTRY ♦ 1 

IF(PBAR.GE.l..AND.NENTRY2.EQ.O) NENTRY « 1 
IFdOUT.GT.l .AND.SNPO.EQ.l) CALL PftT0UT(2»l) 
IF(PBAR.GE.1..AND.NENTRY2.EO.O) NENTRY - NENTRY ♦ 1 
IF C IBR.EQ.O. AND.X2.GT.XLAST) WRITE (6, 2) 

IF( IBR.EQ.O* AND*X2.GT*XLAST) STOP *END OF NOZZLE* 

C 

C BRANCH BACK FOR S.LT.l 

C 

IF( PBAR.LT.i.) GO TO 10 

IF ( NENTRY2*EQ* 0) X22 * XLST 

1 F ( NENTRY2 *EQ* 0 ) Y22 * POUT(NOUT) 

KENTRY2 - NENTRY2 ♦ 1 
TR - T2W/TC 
C 

C CALCULATE SIGHAX FROM EQUATION (58) 

C 

SIGHAX - (82.27 ♦ 75.6124TR - 256 . 889*TR**2 ♦ 95. 928*TR**3 ) 

1 +1.E-3 

C 

C CALCULATE HG1 USING EQUATION (43) 

C 

CALL CHG(T2W,DUM,HG1) 

C 

C CALCULATE LATENT HEAT OF EVAPORATION 

C USING EQUATION (56) 

C 

CAP L * HGl - 4 . 2E ♦$* T2 U ♦ 17 . 1 1 75 365 8E *6 
C 

C CALCULATE CRITICAL DROPLET RADIUS USING EQUATION (59) 

C 

RSTAR - 2.*SIGMAX/(RH0L*RBAR*T2W4AL0G(PBAR) ) 

CALL CCP(T2W,CPfCPl*CPC) 

GAHMA1 - CP1/(CP1 - CAPR/W( 1 ) ) 

MAVG * 0. 

DO 40 I-ltlO 

MAVG - WAVG ♦ 5 IGMA ( I ) *W( I ) 

40 CONTINUE 

MAVG - WAVG/SIG 
WRAT - (Nd)/WAVG>**2 
C 

C CALCULATE Q USING EQUATION (62) 

C 

U - 2 * * ( G AHMA l - l.)*CAPL+(CAPL/(RBAR*T2U) - 0. 5 ) / (RBAR*T2W 
1 * ( GAMHA1 ♦ 1.)) 

C 

C CALCULATE NUCLE AT ION RATE USING EQUATION (61) 

C 

J - (l./U. ♦ Q))*QC*MRAT*SQRT(2.*SIGHAX/(PI*M1**3))*( RHQ**2 
1 /RHOLJ+EXPP (-4.*PI*SIGHAX*RSTAR**2/ (3,*K*T2W) > 

J - J*($IGHAil)/5IG)**2 

IFdOUT.GT.l. AND. SWPO.EQ.l ) CALL PRT0UT(2i2) 

C 

C BRANCH BACK IF J.LT.JM1N 

C 
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1 

T 

r 


SUBROUTINE PARTS 



c 

3 

l 

c 


NUCLEATIQN AND DROPLET GROWTH 



4 

c 





b 



COMMON/COE F CO /COEF (10, 8, 2) 



6 



COHMON/GAMCUT/GAMBAR 



7 



DIMENSION HI 1 f 400) , SL I ( 400 ) * Y I J ( 400) 



6 



COMMON/S IGCQM/SIGMA(10),S IGtPSUPO 



9 



C0MMON/EGCOM/DELTAT,CAPL,NU,PSI,ORDT, 



10 



1 TRiLTILDAtKN, THETA , LAMBDA, Pft,X I , OMEGA , CAP LAM 



11 



REAL NUtLTILDA ? KN» LAMBDA 



12 



COMMON /I NPT/T0,P0,XLASTt PH 1 1,RHC»RN0 



13 



COMMON /SW I TCH/SWEND*SWPO,SwGO 



14 



INTEGER SWEND'SWPO'SNGQ 



15 



COMMON/ 0UT2 /U2 , RH02 » CA P J 2 , T 2 , P 2 , RS T AR2 ,X2 , A2,5 IG2(11), OUMC 



16 



COMMON/CM 0UT3/NENT 



17 



COHMON/CUT 1/ SO ,HQ, MD QT ,RH0STR tTSTAR 



16 



REAL MDOT 



19 



COMMON/CHOUTl/NENTRY t DT,RHOO,Tl,Hl,Ul, USTAR ,URH0,RH0F1 



20 



DIMENSION SIGMSYdO) 



21 



COMMON/AD JUS T/GAMMAltQCt ALPHA, BETA, DSTAR, AS TAR #DT2* 



22 



1 JMIN,DELX,EPS(8),I0uT,0Tl,C¥RATtARAT» JSTOP , JDB ,CP *CP 1 



23 



REAL JMIN 



24 



COMMON/CO NST/TC » K, CAPR , Ml ,M ( 10 ) ,RHOL , R8 Aft , P I , ALPHAC 



25 



REAL K ,N1 


" 

26 



COMMON I , J ,A (400) ,X ( 4 00 ) ,DELN ( 400 ) ,C AP JX ( 400) , R C 400 ) , 



27 



1 MU (4 00) ,TL ( 400 ) , T ( 400 ) ,P (400) ,RSTAR(400) , 



28 



2 Ul 400) ,S (400) , SL (400 ) ,HL(400) ,RH0(40O) ,CAP JY(400> t 


- 

29 



3 F (400) ,G(400) , YS ( 400 ) , DE L Y ( 400 ) ,DELS ( 400 ) , 



30 



4 RPERH1400) ,TS(400) ,OUTMCH( 400) ,CAPP J( 400) 


■ 

31 



REAL MU 


1 

32 



REAL LB AR 



33 



DIMENSION INDT Y P ( 4 00 ) 


■ 

34 

c 





35 


1 

FORMAT ( 1 PROGRAM STOP.#. EXCESSIVE LOOPING ♦+♦♦*') 


= 

36 


2 

FORMAT ( 1H0 , * BEFORE DBL I T * ) 



37 


3 

FORMAT UH , * J,TOLD,T( J + 1 ) , P ( J ♦ 1 ) » , 1 5 , 3E 16 . 8 ) 



38 


4 

FORMAT (* R$TAR( Jd) ,CAPJ-( J),CAPJMJ) • ,3E16.8> 



39 


5 

FORMAT ( 1 DELN( J) ,U(J*1) ,TS( J+l) ,CP1( J+l) S4E16.8) 



40 


6 

FORMAT ( ■ F( J*l) ; C( Jd) ,YS( Jd ) * OELY ( J) S4E16.0) 


- 

41 


7 

FORMAT { 1 HL ( Jd ) , SL ( J d ) , DELS { J ) t S ( Jd ) ■ , 4E 16 .8 ) 



42 


6 

FORMAT ( * RHO ( J d ) « A ( J + 1 ) 1 ,2E16.6) 



43 


9 

FORMAT (’ <R(I, Jd >, I -1,J> *,6816.6) 



44 


10 

FORMAT (' (TL( I, J*l) , I-i, J) S6E16.8) 


= 

45 


11 

FORMAT (» (SL( I, Jd),I«l, J) 1 ,6E16.8) 



46 


12 

FORMAT ( 1 ( HL ( I»Jd)tld»J) ' , 6E 16. 8 ) 


* 

47 


13 

FORMAT ( 1 (Y(I,J+1),I-1,J) , ,6E16.8) 


= 

4b 


14 

FORMATUHO,' AFTER DBLI T* ) 

; 


49 



MAX J - 1000 

- 


50 



JTOP - MAX J 



51 



00 20 I J * 1 , 400 



52 



RPERM(IJ) - 0. 

~ 

■ 

53 


20 

CONTINUE 


“ 

54 

c 



i 


55 

c 


TAKE VALUES FROM PART 2 AS A STARTING POINT 

_ 


56 

c 





57 



T ( 1 ) - T 2 

- 


58 



P(l) - P 2 

' 


59 



RHO(l) - RH02 

: 


60 



U(l) • U2 


3 

61 



SU) - SO 

c 

"■* 

62 



A(l) - A2 

* 


63 



OUTMCHd) * OUMC 

l 

= 

64 



X(l) - X2 



65 



RSTAR( 1) - RSTAR2 



66 



YSd) * 0. 

~ 


67 

c 




— *■ 

68 

c 


CHEMICAL COMPOSITION REMAINS CONSTANT THROUGHOUT 

■ 

' 

69 

c 


PART 3 EXCEPT FOR FORMATION OF LIQUID HATER 

- 


70 

c 


SAVE PART 2 RESULTS 



71 

c 



I 


72 



00 30 IJ-1,10 



73 



SIGMSVdJ) - SIGMA(IJ) 

_ 

— - 

74 


30 

CONTINUE 

: 


75 



SIGSV - SIG 

4 

~ 

76 



J * 1 

- 


77 



I CHG - 0 


"_ 

78 


40 

CONTINUE 

7 


79 



T( Jd) • T ( J) 



80 



P( Jd) - P( Ji 

- 


81 

c 



- 


82 

c 


EQUATION (18) 

- 


46 



original face is 


83 

84 

C 


rvn 

00 50 IJK-2,10 Vrr ' 

POOR 

85 



SIGMA(IJK) - 5 IGMSV ( I JK )/( 1. - Y5(J)) 


86 


50 

CONTINUE 


07 

C 




8b 

C 


EQUATION (17) 


89 

C 




90 



SIGMA(l) « SlGMSV(l) - Y5(J)/0.018 


91 



SIGMA(l) - SIGMA( 1)/(1. - YS(J)) 


92 



SIG - SIGSV - YS(J)/0.016 


93 



SIG • SIG/(1. - Y5 ( J) ) 


94 



CALL SATTEM(TS(J*1),P( J ♦ X ) ) 


95 



U( J*l) - U( J) 


96 

C 




97 

C 


INCREASE THE SIZE OF DELTA X FOR POSITIONS 20 


98 

C 


STEPS PAST THE DROPLET FORMATION CUTOFF 


99 

c 




100 



IF ( ICHG.LT.l 5. AND, ( J-20) .GT.MAXJ) DELX - 1.08*DELX 


101 



IF( ( J-20) .CT.MAXJ) ICHG - ICHG ♦ 1 


102 

c 




103 

c 


INCREASE X BY DELTA X* EQUATION (106) 


104 

c 




105 



X( J+l) - X( J) ♦ DELX 


106 

c 




107 

c 


IF AT END OF NOZZLE* CALCULATION IS FINISHED 


108 

c 




109 



IF(X( J*l> .GT.XLAST) GO TO 240 


HO 



FRCBND - 1. 


111 



IF(J.EQ.l) GO TO 60 


112 



FRCBND - (X(J> - X( J-1))/(X(J*1) - X(J)) 


113 


60 

CONTINUE 


114 

c 




115 

c 


GET THE CRCSS SECTIONAL AREA FOR THE NEW 


116 

c 


X POSITION, EQUATION (107) 


117 

c 




118 



CALL A V 5 X ( XI J+1)*A(J+1) ) 


119 



CALL SATTEM 1 TS ( J ) *P ( J ) ) 


120 



PINF - E XPP (55* 897 - 6641.7/TCJ) - 4. 48 64*ALQG < T (J > )) 


121 



PBAR - (P( J)/PINF)*(SIGHA( 1J/SIG) 


12 2 



TR - T ( J ) / TC 


123 



SIGMAX - (82.27 ♦ 75.612*TR - 256 .889 *TR»*2 ♦ 95. 928 *TR**3 > 

124 


] 

L * 1 # E- 3 


125 



RSTAR(J) - 2.4SIGMAX/(RH0L*RBAR*T( J ) *ALOG( PBAR ) ) 


126 



R(J) - 1 . 001*R5TAR ( J ) 


127 



CALL GRLPRP 


128 



IF(J.NE.l) GO TO 80 


129 



CALL AVSX( XU)-DELX,A JM) 


130 



CAPGAM - MDOT/AJM 


131 



TAUOLD - T(J) 


132 


70 

CONTINUE 


133 



CALL CDEN( RHOF , SO* TAUO LD ) 


134 



CALL CHG(TAUOLD*HGF,DUH) 


135 



CALL CCP(TAU0LD*CP*CP1*CPC) 


136 



UF - SORT ( 2 • * ( HO - HGF ) ) 


137 



FF » RHOF^UF - CAPGAM 


138 



DRHODT - (CP - C AP R’4 S I G ) *RHOF/ TAUOLD 


139 



DUDT - -CP/UF 


140 



FFP - UF4DRH0DT ♦ RHOF* DUDT 


141 



TAUNEW - TAUQLD - FF/FFP 


142 



DEL - ABS(TAljNEW-TAUQLD) / TAUNEW 


143 



TAUOLD - TAUNEW 


144 



IF(DEL.GT.EPS(6) ) GO TO 70 


145 



T JM 1 - TAUNEW 


146 



P JM I * RH0F*CAPR*5IG*TAUNEW 


147 



GO TO 90 


148 


60 

CONTINUE 


149 



TJM1 - T(J-l) 


150 



P JM 1 - PlJ-1) 


151 


90 

CONTINUE 


152 



T JM 1 - (TJM1 ♦ (FRCBND - 1 • ) *T( J) ) /FRCBND 


153 



PJM1 - (PJMI ♦ (FRCBND - 1 • )*P ( J ) )/F RCBND 


154 



T JM - 0 • 5* ( T JM 1 ♦ T ( J ) ) 


155 



P JM « 0.5*(PJM1 ♦ P(J)) 


156 

c 




157 

c 


GET NUCLE AT ION RATE AT J - 1/2 


158 

c 




159 



CALL NUCRAT(CAP JX( J),TJM,P JM) 


160 



DO 100 I-1,J 


161 



RPERMCI > * R ( I ) 


162 


100 

CONTINUE 


163 



NLOOP - 0 


164 


HO 

CONTINUE 


165 



NLOOP - NLOOP ♦ 1 


166 



I F(NLDOP.GT.400) WRITE(6*1) 
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IF(NLOOP.GT.40Q) GO TO 240 
TJP - 0.5*(T(J) ♦ TU*i)) 
PJP - O. 5*( P ( J ) ♦ P( J+l)> 
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16 6 

169 

170 

171 
177 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 
183 
18 4 

185 

186 

187 

188 

189 

190 

191 

192 

193 


196 

197 

198 

199 

200 
201 
202 
203 
2 04 

205 

206 

207 

208 
*09 
210 
211 
212 

213 

214 

215 

216 

217 

218 
219 
2 20 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 
2 37 
236 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 


C 

C GET NUCLEATION RATE AT J ♦ 1/2 

C 

CALL NUCRAT(CAPJY(J)*TJP*PJP) 

DELN(J) - 0. 

CAPPJIJJ - 0. 

IF (CAP JYl J).LT.1») GO TO 120 
IF(CAPJX(J).LT.l.) GO TO 120 
C 

C CALCULATE MEAN NUCLEATION RATE, EQUATION (108) 

C 

CAPPJ(J) - tCAPJYCJ) - CAPJXU) )/ALOG(CAPJY(Jl/CAPJX(J)) 

C 

C CALCULATE NUMBER OP OROPlETS FORMEDt EQUATION (109) 

C 

OELN(J) • CAPPJ( J)*A(J)*DELX/MDOT 
120 CONTINUE 
C 

C SET FLAG TO INDICATE DROPLET FORMATION HAS STOPPED 

C 

IF(DELN (J ) .EQ. 0.. AND, MAXJ.EQ. 1000) MAXJ • J 
C 

C CALCULATE R* FOR DROPLET GROWTH 

C EQUATIONS (57) - (59) 

C 

PINF - EXPP ( 55 • 89 7 - 664 1 * 7/T ( J *1 ) - 4 . 48 64*ALGG ( T ( J ♦ 1 ) ) ) 
PBAR * (P( J + U/PINF )4(SIGMA(1)/5IG) 

TR - T(J*1)/TC 

SIGMAX - (82*27 ♦ 75.612*TR - 256 . 88 9* TR** 2 ♦ 95* 928 *TR**3 ) 
1 41.E-3 

RSTARU+i) * 2.*$IGMAX/(RH0L*RBAR*T( J + 1 ) *ALOG (PBAR ) ) 

C 

C CALCULATE THE TIME 5TEP CORRESPONDING TO THIS 

C DELTA X» EQUATION (110) 

C 

DELTAT - 2.*DELX/(U( J) ♦ U(J*1J) 

DO 130 I«1»J 
R(I) - RP ERM ( I ) 

130 CONTINUE 

CALL GRLPRP 
DO 140 I«ltJ 

I F(NL00P.EQ*1 ) INDTYP(I) * 27 
140 CONTINUE 

C SET DROPLET RADIUS TO ZERO FOR POSITIONS PAST 

C THE DROPLET FORMATION CUTOFF 

C 

IF( J.GT.MAXJ) R { J ) » 0. 

JTOP - J 

IF (J.GT.MAXJ) JTOP - MAXJ 
DO 160 I ■ 1 i JTOP 

IF ( R( I ) /RSTAR( J) . LT. 1 . 1 .ANO .NLOOP. EQ.l ) INDTYP(I) ■ 26 
IF( INDTYP( I ) .EQ.26) GO TO 150 
Z - 1,0001 

C 

C CALCULATE DROPLET GROWTH FOR LARGE 

C DROP LETS t EQUATION (72) 

C 

IF(NLOOP.EQ.l) 

1CALL GRL(Z) 

IF(NLOOP.EC.l) 

1 R ( 1 ) - Z*RSTAR ( J ) 

I F(NLOOP.EQ.l) RPERM(I) - R(i) 

GO TO 160 
150 CONTINUE 

IF(R(I J.EQ.O.) GO TO 160 
Z - R ( I ) /RSTAR l J ) 

C 

C CALCULATE DROPLET GROWTH FOR SMALL 

C DROPLETS, EQUATION (79) 

C 

2 Z - GRS(Z) 

IFiZZ.LT.O.) 11*0. 

R(I) - ZZ*RSTAR( J*ll 
160 CONTINUE 

R ( J 4 I ) - 1,001*RSTAR{ J*l) 

RHQ ( J ♦ 1 ) - P(J*l)/( CAPR+S IG+TIJ+l) ) 

F( J*l) - 0. 

JP1 - JTOP + 1 

CALL NUCRAT (CAPJP1,T(J*1) *P(J^1)) 

DELN(J>1) - A( J41)40ELX4CAPJP1/MD0T 
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260 
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307 
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328 

329 
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334 


ORIGINAL PAGE IS 
c app j ( j ♦ i ) - cap jp i OF POOR QUALITY 

YS4J+1) - 0. 

00 170 I - 1 ♦ J P 1 
C 

C CALCULATE MASS OF DROPLETS FOR TYPE I 

C 

MUU) - 44./3.)*Pl*RH0L*R(l ) **3 
C 

C CALCULATE MASS RATIO OF LIQUID HATER 

C FOR TYPE I, EQUATION <1171 

C 

Y I J ( I ) - MUU)*DELNI I) 

C 

C CALCULATE TOTAL MASS RATIO OF LIQUID 

C HATER, EQUATION (118) 

C 

Y S ( J + 1 ) - Y S 4 J ♦ 1 ) + YIJ(I) 

170 CONTINUE 

DELY(J) * YS ( J* 1 ) - YS ( J ) 

HL4J+1) - 0. 

SL4J+1) - 0, 

JPP1 « J ♦ 1 
DO 200 I-1,JPP1 
TL4I) - T ( J ♦ 1 ) 

IF(I.CT.JPl) GO TO 200 

IF4R4I) .LT.RSTAR4 J+l)> GO TO 190 

KN « L T I LDA / ( 2 * *R ( I ) ) 

DELTA - X I *KN/P R 

CALL CCP(T( J) *CP,CP1,CPC) 

WAVG - 0. 

W C - 0. 

DO 180 IJ-1,10 

WAVG * HAVG ♦ W ( I J) + SI GMA ( I J) 

IF(IJ.NE.l) WC - HC + H ( I J ) *S IGMA 4 1 J ) 

180 CONTINUE 

HAVG - HAVG/SIG 

WC * WC/( SIG - S I GMA 4 1 ) ) 

WRAT - H ( 1 ) /HAVG 

Y 1 - SIGMA(1)/SIG 

GAMMA 1 - CP1/CCP1 - CAPR/HtlJ) 

GAMMAC • CPC/(CPC - CAPR/WCI 

FFF ■ Y1*$GRT( WRAT) ♦ (1, - V 1) *SQRT (HC/wAVG) * 4 GAMMA1/GAMMAC i 
1 *4 (GAMMaC + 1 • ) /( GAMM Al ♦ 1 . >) * 4 CPC/C ?1 1 * ALPHAC 
DELTA • DELTA/ 4XI*KN/PR ♦ 

1 FFF/ l 1 . ♦ 2* *BETA*KN) ) 

TL(I) - T C J ♦ 1 1 ♦ (1./C1. - NU*DELTA) )♦ 

1 ( 1 . - RSTAR4J + l)/RCin*(TS( J*l) - T4J+1)) 

190 CONTINUE 
C 

C CALCULATE ENTHALPY OF LIQUID, EQUATION (1191 

C 

HLI(I) - 4.2E*C3*TLm - 17 ,11 75 3650 E+Q6 
HL4J+1) - HL4J+1) ♦ Y I J C I ) + HL 1(1) 

C 

C CALCULATE ENTROPY OF LIQUID, EQUATION (120) 

C 

SLI4I) - 4.2E+03*AL0G(TL( I) ) - 2 . 001 4962 81E *04 
SL ( J ♦ 1 ) - SL ( J + 1 ) ♦ YIJ(l)*SLl(I) 

200 CONTINUE 

TJ - 0,5*(TIJ*1) ♦ TtJ)) 

PJ - 0,5*4 P( J+l ) ♦ P(J)) 

CALL 5ATTEM4TS J,PJ> 

JZ • 1 

IF4TSJ.GT.100d. > JZ - 2 
CALL CHG ( TS J , DLM , HGX) 

L BAR - HGX - 4 »2E+3*T S J ♦ 17. 1175 3658E *6 
JZ - 1 

IF4TJ.GE.100C.) JZ - 2 
CALL CCPlTJfDUMfCPZtDUM) 

C 

C CALCULATE CHANGE IN ENTROPY, EQUATION (83) 

C 

0ELS4J) - CLBAR - CPZ*4TSJ - TJ))*DELY(J) 

1 *(l./TJ - l./TSJ) 

C 

C CALCULATE ENTROPY OF THE MIXTURE, EQUATION (124) 

C 

S 4 J ♦ 1 ) - S( J) + DELS 4 J ) 

TOLD - T ( J* 1 ) 

IF4J.LT. JOB) GO TO 210 
HRITE46,2) 

HRITE(6,3> J,T0LD,T4 J + l) , P 4 J + 1 ) 

WRITE ( 6, 4 ) RSTAR4 J+l) , CAP JX ( J ) ,CAPJY(J ) 

WRITE (fa, 5) DELN4 J),U( J + l) ,TS( J + l) ,CP1 
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335 

336 

337 
336 
339 
3*0 
3*1 
3*2 
3*3 
3** 
3*5 
3*6 
3*7 
3*8 
3*9 

350 

351 

352 

353 
35* 

355 

356 

357 

358 

359 

360 

361 

362 

363 
36* 

365 

366 

367 

368 

369 

370 

371 

372 

373 
37* 

375 

376 

377 

378 

379 

380 

381 

382 

383 
38* 

385 

386 

387 
386 

389 

390 

391 

392 

393 
39* 

395 

396 


WRITE* 6*6) Ft J*l),GCJ+l) , YS ( J ♦ 1 ) *DEL Y t J ) 

WRITE* 6*7) HL(J4l)*SL( J*l> *DEL$U) *S(J*1> 

W RI TE I 6 ,8 ) RhO C J ♦ 1 1 * A t J ♦ 1 ) 

WRiTE(6*9) (R(I JK ) , I JK-1, J ) 

WRITE *6, 10 I (TLUJK)*IJK-1*J) 

WRITEt6*ll> (SLI (I JK),I JK-1, J I 
WRITE (6*12) IHLltl JK), I JK«1*J) 

WRITE (6,13) (YIJ< I JK) ,1 JK-1, J> 

210 CONTINUE 
CALL DBLIT 

CALL CCP ( T ( J* 1 ) *CPS*CP1*CPC) 

GAMMAS - CP1/(CP1 - CAPR/WCD) 

OUTMCH ( J*1 ) - GAMMAS*CAPR+SIG*Tt J ♦ 1 1 
OUTMCH(JU) • U( J>1)/SQRT(0UTMCH( J4l) ) 

IF t J *LT, JD B ) GO TO 220 
WRITE(6,1*» 

WRITE ( 6* 3 ) J*TQLD*T< J + l) *PU+1) 

WRITE *6**) RSTARt J ♦ 1 ) ,CAPJX(J) *CAPJY(J) 

W'RITE(6*5) OELNt J) *U( J*l) ,TSl J+l) *CP1 
WRITE *6*6) Ft J*l> *G* J*l> *YSC J ♦ 13 *DELY t J ) 

WR ITE t 6* 7 ) HL*JU)*$L(J>1) *DELS( J) *S(J*1) 

W R I TE * 6 , 8 ) RHO ( J + l ) , A t J +1 ) 

WRITE*6*9) *RtI JK) *1 JK-1* J) 

WRITE *6*10) tTLtl JK) *1 JK-1, J) 

WRITE (6,11) t SLI (I JK) * I JK-1* J ) 

WR ITE * 6* 12 ) *HLI tIJK) * I JK-1 * J ) 

WRITE *6,13) t Y I Jt IJK), I JK-1*J) 

220 CONTINUE 
C 

C CONVERGENCE TEST 

C 

DEL - (T(J*1> - TO LD ) / T OLD 
TOLD - T ( J + l I 
C 

C TEST FOR TEMPERATURE CONVERGENCE EQUATION *133) 

C IF NO CONVERGENCE GO BACK TO DROPLET GROWTH 

C 

IF(ABStDEL) *GT*EPS(7) ) GU TO X 10 
C 

C CONVERGENCE, » .OUTPUT RESULT5 AND GO ON 

C TO THE NEXT VALUE OF J 

C 

WRITE* 7) J,*RU>tMum,HLim*5lIU),TLm*YIJ*I),I-l.J) 
GAMBAR - AIOG<T(J*1>/TO)/ALOGCPIJ*XI/PO> 

GAMBAR - 1.7(1. - GAMBAR) 

IF( I0UT.EQ,2»AND.SNP0.EQ.l ) CALL PRT0UT(3*1) 

UO 230 I - 1 * J 

IF(INDTYPtl) »EQ« 26 ) RP ERM C I ) - R 1 1 ) 

230 CONTINUE 

GAMBAR - ALOGt T(J+1)/TQ)/AL0G* PtJ+l)/PG> 

GAMBAR - l./tl. - GAMBAR) 

WRITE* 33) J * X ( J ) *DELN ( J ) , t RtKK ) *KK-1, J) 

J - J ♦ l 
C 

c IF AT TERMINAL J* CALCULATION IS FINISHED 

C 

IFt J.GT.JSTOP) GO TO 2*0 
GO TO *0 
2*0 CONTINUE 

IF(SWPO.EO.l) CALL PRT0UT(3*0) 

RETURN 

END 
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SUBROUTINE PLTQUT 


1 

2 

3 

4 
b 
6 

7 

8 
9 

10 

11 

12 

13 

14 

15 
lb 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 
27 
26 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 


C 

C SUBROUTINE PLTQUT GENERATES PLOT OUTPUT OF PROGRAM 

C RESULTS. SEVERAL SUPPORTING SUBROUTINES ARE A PART 

C OF THE LANGLEV GRAPHICS SYSTEM. THESE INCLUDE* 

C ASCALE 

C AXES 

C CALPLT 

C CHARST3 

C CHARST4 

C CHARS 12 

C DASHPLT 

C LINPLT 

C NFRAME 

C NOTATE 

C PNTPLT 

C REPCHAR 

C FOR MORE INFORMATION ON THESE ROUTINES SEE M L ANGLEY 

C GRAPHICS SYSTEM”* CENTRAL SCIENTIFIC COMPUTING COMPLEX 

C DOCUMENT G-3. 

C 

COMMON I , J,A(400) *X(400) *DELN(400) *CAP JXtAOO) *R<400) * 

1 MU (400) *TL(400) *T(400) * P ( 400 ) . RSTAR ( 400 ) * 

2 U ( 400 ) *S (400) *SL(400) ,HL(400) *RH0(400) *CAPJY(400) t 

3 F(400) *G(400) ,YS(400) » DE LY ( 400 )* DELS ( 400 ) * 

4 RPERM(4C0) * T S ( 400 ) , OUTMCH ( 400 ) , CAPP J ( 400) 

REAL HU 

COMMON/ OUT 1/SO* HO *MDOT * RHOS TR * TS TAR 
REAL MOOT 

DIMENSION OUM ( 400 ) 

COMMON/ 5VX2/X22*Y22 

C0MHON/EXPDAT/NEXP,XEXP(50) *PEXP (50) ,PREF 
COMMON/P LTB LK/XM IN, XM AX* PM IN ,P M AX , TM IN , TMAX , 

1 YMIN, YMAX*DTMIN,DTMAX,DNMIN,DNMAX,RMIN,RMAX, 

2 PRMIN,PRMAX,TLMI N, TLMAX , RHOM IN ,RHOMAX , 

3 SMIN*SMAX*UMIN*UMAX*XMMIN,XMMAX,!PSLCT( 15) 

COMMON/ INPT/TO*PO*XLAST*PHI1*RHC,RNO 
COMMON/PLTOUT/NOUT*TOUT(bOO) , POUT (600 ) * XOUT (600 ) , 

1 DELTT ( 600 ) *UOUT(600) *0UTMAC(600) 

DIMENSION YP ( 400 ) 


41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 


01 ME NS ION P A T ( 2 ) 

DIMENSION PATX (4 ) ,PATY (4 ) 

DATA PAT/0.5,-0.25/ 

DATA PATX/O. tl.25,0.,1./ 

DATA PATY/O. 58,0. 58,0. ,1./ 

C 

CALL CHARS12 
CALL REPCHAR(2,12 J 
CALL CHARST4 
CALL CHARST3 
NP - J - 1 
XSV - XII) 

X<1) - 0. 

CALL ASCALE(X,8.,NP,1,10.) 

IF(XMIN,NE.O..OR.XMAX.NE.O.) X(NPU) - XMIN 
1 F C XMIN.NE.O. .OR.XMAX.NE.O. ) X(NP*2) - (XMAX - XMINI/8. 
XOUT (NOUT+1 ) - X ( NP ♦ 1 ) 

XOUT ( NOUT +2 ) - X ( NP +2 ) 

XU) - XSV 

IF(IPSLCT(2) • EQ .0 ) GO TO 40 

C 

C PRESSURE RATIO VS X 

C 

NOUTM - NOUT - 1 
DO 30 IJ-1,NP 
DO 10 I JK ■ 1 , NOUTM 

1F(X(IJ>.GE.X0UT(IJK).AND.X(1J).lE.XUUT(IJKU)) GO TO 20 
10 CONTINUE 
I JK - NOUTM 
20 CONTINUE 

YP(IJ) - POUT(IJK) ♦ (XUJ) - XOUT(IJK)) 

1 * ( POUT ( I JK ♦ 1 ) - POUT( I JK) )/(XOUTU JK ♦ i) - XOUT(IJK)) 
YP(IJ) - P(IJ)/1 YP(IJ)*PO) 

30 CONTINUE 

CALL ASCALE(YP,6.,NP,l,10.i 

IF (PRMIN.NE.O.OR.PRMAX.NE.O) YP(NPU) - PRMIN 
IF(PRHIN.NE.O.OR.PRMAX.NE.O) YP(NP*2) - ( PRMAX - PRMINJ/6. 
CALL AXES(0.tO««0.»B. 9 XCNP*l)»X(NP+2l *1.,0. • 

1 1HX, 0. 2 ,-l ) 

CALL AXES (0. ,0*,90. ,6* *YP (NP+1) ,YP ( NP + 2) ,l.,0., 

1 7H P RATIO, 0.2, 7) 

CALL LINPLT(X,YP,NP, 1,0, 0,0,0) 
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03 

0*. 

05 

66 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 
15B 

159 

160 
161 
162 

163 

164 

165 

166 


C 

c 

c 


XP LOT - ( X2 2 - XiNPnn/XlNP*21 
YPLOT - 0, 

CALL PNTPLTCXPL0T,YPL0T.4,1> 
CALL NF RAM E 

CALL CALPLT C4.5tl.5t -3) 

40 IF l IP5LCT (7) * EU • 0 ) GO TO 60 

CRITICAL RADIUS V5 X 


DO 50 U-l.NP 
YPUJI * RSTAR(IJ) 

50 CONTINUE 

CALL ASCALEC YP t 6. t NP tl 1 10 • ) 

IFIRMIN. NE.O.. OR. RMAX. NE.O.) YP(NP*1) - 
IF CRHIN.NE. 0. .GR.RMAX.NE. 0 . ) YPCNP+2) • 
CALL AXES(0.»0.,0.»8.,X(NP + l)*XCNP«>2).l 


RHIN 

( RM AX-RM I N ) /6 « 

• t 0 • t 


1 1HX » 0. 2 .-1 ) 

CALL AXES ( 0. .0 . .90. »6. 


YP ( HP+l) .YPCNP+2 ) *1. *0. » 


1 5HRSTAR.0.2.5) 

CALL L INPLTCX. YPtNP.l.O.O.OtO) 


CALL NFRAME 

CALL CALPLTC4.5fl.5t-3> 

60 IFUPSLCTC8).EO.O> GO TO 140 


C 

c DROPLET RADIUS (AND R*J V5 X 

C 

N RAD * NP/10 

NRAD - NR AD - 1 

IF C N RAD* GT • 10 ) NRAD * 10 

IRAD - 5 

DO 120 JRAD* 1 » NRAD 
REWIND 7 

RE ADC 7 ♦ E ND* 80 ) JJJ»(DUM(IJ).DUMl.DUM2»DUM3iDUM4tDUM5?IJ*ltJJJ> 
80 IFCEQF (7) .NE.O) GO TO 90 
IFC JJJ.LT. IRAD) GO TO 70 
YPCJJJ) - DUMCIRAD) 

GO TO 70 
90 CONTINUE 

NPP * NP - IRAD ♦ 2 
YPCIRAD-1) - 0. 

DO 100 IJ-l.NPP 
KRAD - IRAD ♦ IJ - 2 

ir-(YP(KRAD).LE.10.**RMIN) YP(KRAD) - 10.**RMIN 
YP(KRAD) - ALOGIOC YP CKRAD ) ) 

100 " CONTI NUE _ m 

IFURAO.EO.il CALL ASCALEC YPCIRAD-1) f 6. tNPP tit 10. ) 

IFC JRAD.NE.l) GO TO 110 

IFC RHIN. NE.O. .OR. RMAX, NE.O.) YPCNP+l) - RMIN 
IFCRMIN. NE.O. .OR. RMAX. NE.O.) YPCNP+2) - C RMAX-RH IN I /6. 

110 CONTINUE 

IF C J RAO. EO. 1 ) 

1C ALL AXES(0.»0.»0.»8,.X(NP*l)tX(NP*2)tl.*0.t 

1 • D C ISTANCE ) A ( LONG ) N{ OZZLE X* M)».0.2t-32l 
IFC JRAD. EO. 11 

ICALL AXES (0. tO • .90. .6 . » YP C NPU) t YP (NP + 2) 1 1 .5. 0. t 

2 1 DC ROPLET )RCAD1US» I2LS30G RSUISN AND S2LS30G R*D*SNi H . 

3 0.2.55) 

NPP - NPP - 1 

IFCYPC IRAD) .NE. RHIN) Y P < IRAD) - ALOGIOC RSTARC IRAD) ) 

CALL LINPLTCXCIRAD) .YP(IRAO) tNPPflfOtOfOtO) 

IRAD - IRAD ♦ 10 
120 CONTINUE 

DO 130 IJ-l.NP 
YPUJ) - RSTARC I J) 

IFCYPCI J).LE.10.**RHIN) YPUJ) - 10.**RHIN 
YPCIJ) - AL0G10CYPC I J) ) 

130 CONTINUE 

CALL LINPLTCXfYPtNPf lfOtOiOtO) 

CALL NFRAME 

CALL CALPLTC4.5fl.5t-3) 

140 IFCIPSLCTC9) .EQ.O) GO TO 280 
. C 

c TEMPERATURE OF LIQUID VS X 

C 

DO 150 IJ-l.NP 
YPCIJ) - TS C I J I 
150 CONTINUE 

NRAD - NP/10 
NRAD - NRAD - 1 
I F ( NRAD . GT . 10 ) NRAD - 10 
IRAD - 5 

DO 200 J RAD" 1 * NR AD 
REWIND 7 
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167 160 CONTINUE 

1 6B READ* 7, END- 170) J J J 9 ( DUM4 , UUM 1 , DUM2 t DUM3 t DUM ( I J ) ,DUM5 . 1 J-l * J J J > 

16? 170 IF ( EOF C 7 ) .NE.O) GO TO 180 

170 IF* JJJ.LT.IRAD) GO TO 160 

171 YP* J J J ) * DUM ( IRAQ) 

172 GO TO 160 

173 180 CONTINUE 

174 NPP - NP - IRAD ♦ 2 

175 YP(IRAD-l) - 0. 

1^6 IF ( JRAD.EQ.l) CALL AS C ALE < YP (I R AD-l ) * 6 . » NP P ,1 f 10 . ) 

177 IF(JRAD.NG«1I GO TO 190 

178 IF (TLMIN.NE.C. .OR.TLMAX.NE.O. ) YP*NP*1| - TLMIN 

179 IFtTLMIN.NE .0. .OR.TLMAX.NE.O. I YP(NP*2> - { TLMAX - TLMIN>/6. 

180 190 CONTINUE 

181 IF{ JRAD.EQ.l) 

182 1C ALL AXE5(0.'0.tO.»6.*X(NP4>l)«X(NP+2) . 

183 1 1HX * 0. 2 t -1 ) 

184 IF (JRAD.EQ.l ) 

18 5 1C ALL AXES (0. tO .»90. *6. *YP (NP*1) «YP (NP«2) * I., 0. 9 

186 1 1HT, 0.2*1) 

107 NPP - NPP - 1 

188 CALL L INPLT*X ( IRAD ) ♦ YP ( I R AD ) »NPP 1 1 , 0 > 0 ,0 , 0 ) 

189 IRAD - IRAD ♦ 10 

190 200 CONTINUE 

191 DO 210 IJ«1,NP 

192 YPUJI * T ( I J ) 

193 210 CONTINUE 

194 CALL LINPLT(X,YP,NP,lt0t0,0»0) 

195 DO 220 IJ-l f NOUT 

196 TOUT* I J ) - T OU T * I J ) *T0 

197 220 CONTINUE 

198 CMAX - YP ( NP ♦ 1 ) ♦ 6.*YP(NP*2) 

199 DO 230 IJ-IyNOUT 

200 IF < TOUT ( I J ) .LE.CMAX) GO TO 240 

201 230 CONTINUE 

202 IJ - 1 

203 240 I8EG - IJ 

204 DO 250 IJ-l.NOUT 

205 IF ( X OUT ( I J) • GE • 20. ) GO TO 260 

206 250 CONTINUE 

207 IJ - NOUT 

208 260 IEND « IJ 

209 NOUTX • IEND - I&EG > 1 

210 TSAV1 - T CUT ( I END* 1 ) 

211 TSAV2 - TOUT ( X END+2 ) 

212 T OUT * I END+ 1 > « YPlNP+1) 

213 TOUT t I END + 2 ) • YP(NP*2I 

214 TM>5 - XOUT ( 1 ENO + 1 ) 

215 TMP6 - XOUT C I E ND ♦ 2 ) 

216 XOUT ( I END+1 J - XQUTtNOUT+l) 

217 XOUT ( I END + 2 ) • XOUT ( NOUT+2 ) 

218 CALL DASHPLTtXOUT (1BEG) »TOUT(IBEC) *NQU TX 1 1 * PAT »2 ) 

219 XOUT ( IEND+1 ) - TMP5 

220 XOUT ( I END *2 ) - TMP6 

221 TOUT ( I END+1 ) - TSAV1 

222 TOUT ( I ENO+2 ) - TSAY2 

223 DO 270 IJ-1,N0UT 

224 TOUT* I J) - TOUT * I J ) / TO 

225 270 CONTINUE 

226 CALL NFRAME 

227 CALL CALPLTI4.5tl.59-3) 

228 280 IF* I PSLCT ( 10) • EQ . C ) GO TO 300 

229 C 

230 C DENSITY VS X 

231 C 

232 DO 290 IJ-l.NP 

233 YP(IJ) - RHO(IJ) 

234 290 CONTINUE 

235 CALL ASCALE(YP,6.»NP,1,10.) 

23b IFiRHOMIN. NE.O.. OR. RHOMAX. NE.O.) YP<NP*1) - RHOMIN 

237 IF(RHOMIN. NE.O.. OR. RHOMAX. NE.O.) YP(NP*2) * (R HON AX - RH0MIN)/6. 

238 CALL AXES(0. *0 . 9 0 . * 8. , X t NP ♦ 1 ) * X * NP*2 > 1 1 . . 0 . , 

239 1 lHXf0.2»“l) 

2 AO CALL AXES *0. *0. , 90 . . 6 . , YP * NP ♦! ) , Y P ( NP+2 ) tl. ,0. » 

241 1 3HRH0 1 0 • 2 9 3 ) 

242 CALL LIKPLTUtYPtNPtliOtOfOtO) 

243 CALL NFRAME 

244 CALL CALPLT(4.5,1.5t-3) 

245 300 IF* IPSLCT*13).EQ.O) GO TO 340 

246 C 

247 C MACH NUMBER VS X 

248 C 

249 CALL ASCALEtOUTMACf6.9NOUT.l9lO.) 

250 lFtXMMIN.NE.O. .OR.XHMAX.NE.O. ) OUTM A C * NOUT* 1 ) -XMNIN 
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251 IF(XMMIN.NE.O. .OR .XMM AX .NE . C. ) QUTMAC ( NQUT + 2) - (XMHAX - XMMlN)/6. 

252 CALL AXES(0. ,0 . ,0 . * 8. ,X< NP* 1 ) , X1NP+2) ,1.,0. , 

253 1 lHX,0.2f-U 

2 54 CALL AXE$(0.,0.,90. ,6 . , OUTMAC < NOUT* 1 > , OUTMAC (N0UT*2 ) , l.,Q., 

255 1 4HHACH,0.2,4) 

256 DD 310 IJ-WNQUT 

257 IFtOUTMAC ( I J ) . GE .OUTMAC ( NOUT* 1) ) GO TO 320 

258 310 CONTINUE 

259 IJ - 1 

260 320 NOUTX - NOUT - IJ ♦ 1 

261 CALL LIhPLT(XQUT(I J),OUTMAC(IJ> ,NQUTX, 1,0, 0,0,0) 

262 DO 330 IJ*1,NP 

263 YPUJ) - OUTMCHdJI 

264 330 CONTINUE 

265 YPINP + l) - OUTMAC i NOUT *1 ) 

266 Y P ( NP *2 ) » OUTMAC ( NOUT *2 ) 

267 CALL LINPLT(X,YP,NP,1,0,0»0,0I 

266 CALL NFRAHE 

269 CALL CALPLT(4.5,1.5,-3) 

270 3 AO IFUPSLCT(l) .EQ.O. AND. I PSLCT(141 .EQ.O) GO TO 400 

271 C 

272 C NORMALIZED PRESSURE VS X 

273 C 

274 00 350 IJ-1»NP 

275 YP(IJ) - P(IJ)/P0 

276 350 CONTINUE 

277 THP1 • YP(NP) 

278 YP(NP) - AM1NKP ( NP )/ PO , POUT (NOUT )) 

279 CALL ASCALE(YP,6.,NP,1,10.) 

280 1F(PHIN.NE.C..0R.PMAX.NE.0.) YP(NP*1> - PMIN 

281 IF (PMIN.NE.O* .QR.PMAX.NE.O. ) YPtNP+2) - (PM AX-PM IN ) /6. 

282 YP(NP) - TMP1 

283 CALL AXES(0. ,0,, 0. ,8. ,X(NP*1) ,X(NP + 2i ,1. ,0. , 

2 8 A 1 1 D ( I STANCE ) A< LONG )N(OZZLE X, M)', 0.2, -321 

205 CALL AXES(0.,0.,90.»6.,YP(NP*1),YP(NP+2),1.,0., 

286 1 , N ( ORMAL I ZED ) S ( TAT IC )P(RESSURE, P/PSDO) 1 * 

287 2 0.2,40 

288 CALL LINPLT(X,YP,NP, 1,0*0, 0,0) 

289 TMP1 - X ( NP 4 1 ) 

290 TMP2 - X ( NP *2 ) 

291 TMP3 - YP(NP*1) 

292 TMP4 - YPINP+2) 

293 XOUT ( NOUT + 1 ) - TMP1 

294 X0UT(N0UT*2) - TMP2 

295 POUT ( NOUT + 1 ) • TMP 3 

296 POUT ( NQUT*2 ) - TMP4 

297 CM AX - YP(NPU) ♦ 6.*YP(NP*2) 

298 DO 360 IJ-1,N0UT 

299 2 F < POUT C I J> .LE.CMAX) GO TO 370 

300 360 CONTINUE 

301 IJ - 1 

302 370 NOUTX - NOUT - IJ ♦ 1 

303 CALL DA$HPLT<XOUT(I J),POUT(I J) , NOUTX , 1 • PAT , 2 1 

304 IF(NEXP.EO.C) GO TO 390 

305 00 380 IJ-1,NEXP 

306 XPLOT - XEXP(I J)/39.37 

307 XPLOT - (XPLOT - X ( NP *1 ) ) /X ( NP *2 ) 

308 YPLOT - PEXP( I JJ/PREF 

309 YPLOT - (YPLOT - YP (NP *1 ) >/YP INP*Z) 

310 CALL PNTPLT (XPLOT, YPLOT, 1, l ) 

311 380 CONTINUE 

312 390 CONTINUE 

313 XPLOT - (X22 - X(NP*1))/X(NP*2 ) 

314 YPLOT - (Y22 - YPCNPU) l/YP(NP*2) 

315 CALL PNTPLTIXPLOT, YPLOT, 901,1) 

316 IF ( I P5LCT (14}.NE.O) GO TO 400 

317 CALL NFRAME 

318 CALL CALPLT (4.5, 1.5, -3) 

319 400 IF(IPSLCT(31. EQ.O. AND. 1PSLCTC14). EQ.O) GO TO 440 

320 C 

321 C NORMALIZED TEMPERATURE VS X 

322 C 

323 DO 410 IJ-1,NP 

324 YPIIJ) - T(IJ)/TO 

325 410 CONTINUE 

326 TMP 1 - YP(NP) 

327 YP(NP) - AMIN1(T(NP)/T0, TOUT (NOUT) I 

328 CALL ASCALE(YP,6.,NP,1,10.) 

329 IFtTMIN.NE.O. .OR.TMAX.NE.O. ) YP(NP*i) • TM IN 

330 IFUMIN.NE.O.. OR.TMAX.NE.O.) YP(NP*2) - (TMAX-TMIN )/6* 

331 YP(NP) - TMP 1 

332 IF(IPSLCT(3).NE.O) 

333 1C ALL AXESfO. ,0.,0.,8. , X ( NP+ 1 ) , X( NP+2 ) t 1. ,0 . , 

334 1 1HX , 0 • 2 , -1 ) 
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335 IF (IP5LCTC3) .NE.Q) 

3 36 1C ALL AXESfO. ,0.,90. ,6. ,YP(NP + 1) ,YP(NP*2) ,1.,0. , 

337 1 AHT/TO *0.2 , A ) 

338 IFUPSLCT(IA).NE.O) 

33 9 1 CALL AXES 1-1. #0. f 90 • »6. • Y P ( NP*1 > f YP! NP*2) ilo 0. » 

3 AO 1 AHT/TO, 0.2, A) 

3 A 1 CALL LIKPLT(X,YP,NP, 1,0, 0,0,0) 

3A2 THP 1 - X { NP ♦ 1 ) 

3 A 3 TMP2 - X ( NP ♦ 2 ) 

3 A A TMP3 - YP ( NP ♦ 1 ) 

3 A 5 TMP A - YP ( N P + 2 ) 

3A6 XOUT ( NOUT *1 ) - TMP1 

3 A 7 XOUT ( NOUT ♦ 2 ) - TMP2 

3 A 8 T OUT ( NOUT + 1 ) - TMP3 

3 A 9 TOUT l NOUT+2 ) - ThPA 

350 CMAX - YP ( NP ♦ 1 ) ♦ 6.*YP(NP*2) 

351 DO A20 I J**l ,NDUT 

352 IFITOUT(IJ) «LE .CMAX ) GO TO A30 

353 A20 CONTINUE 

35A IJ « 1 

355 A30 NOUT X • NOUT - IJ ♦ 1 

35b CALL LINPLT(XOUTIIJ) .TOUT (I J) *NOUTX, 1,0, 0,0,0) 

357 IFI IPSLCTt 1A) • NE . 0 ) CO TO A AO 

358 CALL N FRA ME 

359 CALL CALPLT (A.5,1.5,-3 ) 

360 AAO IFdPSLCT(12).EQ.0.AND.IPSLCT(14) .EQ.O) CO TO 510 

361 C 

362 C FLOW VELOCITY VS X 

363 C 

36A DO A50 IJ-1,NP 

365 YP(IJ) - U(IJ> 

366 A50 CONTINUE 

367 CALL ASCALEt YP , 6 . , NP , 1 ,10 . ) 

368 IFtUMIN.NE. 0* .OR.UMAX.NE.O* ) YPtNP+l) - UMIN 

369 IF ( UMIN.NE. 0. .OR, UMAX. NE .0. ) YP(NP+2> - (UMAX - UMIN)/6. 

370 UOUT (NOUT ♦ 1 ) * YPINP+1) 

371 UOUT { N0UT*2 ) - YP(NP*2) 

372 IF(IP5LCT(12) .NE.C) 

373 1C ALL AXES(0. ,0 . , 0 . , 8 . , X ( NP ♦ 1 ) , X ( NP *2 > tl.,0. » 

37A 1 1HX ,0 • 2 , — 1 ) 

375 IF(IPSLCT(12).NE.O) 

376 1 CAL L AXES (0. ,0.,90. ,6. »YP (NP*1) *YP (NP+2) tl.*0.» 

377 1 1HU, 0.2,1) 

378 IF(IPSLCTUA).NE.O) 

379 1C ALL AXES (-2.,0.,90. ,6. iYP ( N P ♦ 1 1 ,YP (NP*2) ,1 *« 0. , 

380 1 1HU, 0.2,1) 

381 CALL LINPLT (X,YP,NP, 1,0, 0,0,0) 

382 DO A60 IJ*1,N0UT 

383 IF(UGUT( IJ) «GE. (YP(NP+l) + 6.*YP(NP + 2i )) GO T0A70 

38 A A60 CONTINUE 

385 NOUTY - NOUT 

386 GO TO A80 

387 A70 NOUTY - IJ - 1 

388 ABO CONTINUE 

369 DO A90 IJ«1,N0UT 

390 IF(U0UT(IJ).GE.YP(NP+1) ) GO TO 500 

391 A90 CONTINUE 

392 IJ - 1 

393 500 NOUT X - NOUTY - IJ ♦ 1 

39 A XOUT ( NOUT Y* 1 ) - X0UTIN0UT + 1) 

395 XOUT ( NOUT Y+2 ) - X0UT(N0UT*2) 

396 UOUT (NOUTY ♦ 1 ) - YP(NPd) 

397 UOUT ( NOU TY ♦ 2 ) » YP(NP*2) 

39tt CALL LINPLT(XOUT(IJ),UOUTUJ> ,NOUTX, 1,0 ,0, 0 ,0) 

399 CALL NFRAME 

AOO CALL CALPLT (A. 5, 1.5, -3) 

AO 1 510 IF(IPSLCT(4). EQ.O. AND. IP$LCT(15J.EQ.G> GO TO 530 

A02 C 

A 03 C MASS FRACTION VS X 

AO A C 

AO 5 NPM * NP - 1 

A06 DO 520 IJ-1,NP 

A 07 YP( I J) - YS ( I J) 

A 08 520 CONTINUE 

A09 CALL ASCALEIYP ,6. , N P ,1,10. ) 

AiO IFIYMIN.NE.O..OR.YMAX.NE.O.) YP(NPU) - YM I N 

All IF ( YMIN.NE.O. .OR.YMAX.NE.O. ) YP(NP+2) » ( YHAX-YM IN ) ft . 

A 1 2 CALL AXES(0.,0.,0.,8.,X(NP+l),X(NP + 2)»l.,0.t 

A13 2 ■DUSTANCE ) A ( LONG >N(OZZLE X, M)*,0.2,-32) 

A 1 A CALL AXE5(0.,0. ,90. , 6 . , YP ( NP + 1 ) , YP ( N P+2 ) ,1.,0. , 

A15 1 * M ( AS S ) F ( RAC T I ON OF ILdQUIO ) H ( ATER , W) 1 , 

A 16 2 0.2, AG) 

A 17 CALL LINPLT(X,YP,NP, 1,0, 0*0,0) 

A18 IF(IPSLCT(15).NE.O) GO TO 530 


55 



<i 19 
A20 
A21 
A22 
4i 2 3 
A2A 
A 25 
A 26 
A 27 
A26 
A29 
A 30 
A 3 1 
A32 
A33 
A3A 
A 35 
A 36 
A 37 
A 38 
A 39 
A AO 
A A 1 
A A 2 
A A 3 
AAA 
A A5 
A A6 
A A7 
A A 8 
A A9 
A 50 
A 5 1 
A52 
A53 
A5A 
A55 
A56 
A57 
A5B 
A59 
A60 


CALL NFRAME 

CALL CALPLTIA. 5,1. 5,-3) 

530 IFIIPSLCT(5).EQ*Q.AND.IPSLCT(15).EQ.O); GO TO 580 

C 

C TEMPERATURE DIFFERENCE VS X 

C 

DO 5A0 IJ-i.NP 

YP(I J) * TS ( IJ) - T(IJ) 

5 AO CONTINUE 

CALI A$CALE(YP,6.,NP,1,1Q.) 

IF(DTHIN.NE.O..OR.QTMAX.NE.O.) YP(NP*1) - OTMIN 
IFIDTMIN.NE.O. . DR.DTMAX.NE. 0. ) YP(NP^Z) 

1 - { DTMAX - D T Pi IN ) /6. 

IF(IPSLCT(5)*NE*0) 

1CALL AXES (0.*0. *0. >6 .*XINP*1) *X (NP*Z> *1.*0. , 

1 1HX * 0* 2 * -1 ) 

IF ( I PSLCT ( 5) .N E • 0 ) 

1C ALL AXES (0* *0* *90 » *6. *YP(NP*1. J *YPINP*2) * l • * 0 • * 

1 6H0ELTAT * 0 • 2 »6 ) 

IF(IPSLCTtl5).NE.O) 

1CALL AXES (-l.,0.,90. , 6 . , YP l NP+1 . ) , YP l NP *2) *1.* 0. » 

2 ' T ( EMPER ATURE ) D ( I F FE RENCE * ) T$D ( S ) *N ( % ( PSDISNS )- ITIOGSN, K*, 

3 0.2*58) 

CALL LINPLTIX, YP ,NP » 1 * 0 *0 .0 * 0 ) 

DO 550 I J*1 * NOUT 
IF(XOUTUJ).GE.Xm) GO TO 560 
550 CONTINUE 

NOUTX - NOUT 
GO TO 570 
560 NOUTX - IJ 
570 CONTINUE 

DELTT(N0UTX*1) - YPINP+l) 

DELTT ( NOUTX+2 ) - YP(NP*2) 

X0UT<N0UTX*1)-X<NP*1) 

X0UT(N0UTX^2)-XINP^2) 

CALL L I NPLTUOUT.DELTT, NOUTX *1,0, 0,0,0) 

IF(IPSLCT(15).NE.O) GO TO 580 
CALL NFRAME 

CALL CALPLT (A.5*1.5*-3) 

580 IF( IPSLCT (6) .EQ.O.AND. IPSLCT115J.EQ.0) GO TO 610 
C 

C NUCLE AT ION RATE VS X 


A61 
A62 
A63 
A6A 
A65 
A66 
A67 
A68 
A 6 9 
a70 
A71 
A72 
A73 
A7A 
A75 
A76 
A77 
A 78 
A79 
A80 
A81 
A 8 2 
A 83 
A 6A 
A85 
A86 
A87 
A6B 
A89 
A90 
A91 
A92 
A93 
A9A 
A 95 
A96 
A97 
A98 
A99 

500 

501 

502 


DO 590 I J* 1 * NPM 
YP(IJ) - CAPPJIIJ) 

590 CONTINUE 

CALL A5CALE( YP,6. *NPM*1, 10. ) 

IF1DNMIN.NE.0. .OR.ONMAX.NE.O. ) YP(NP) - DNMIN 
IF { DNMIN.NE • 0. .OR .ONH AX.NE .0, ) YP(NP + 1) 

1 - { DNMAX - DNM IN } /6. 

DO 600 IJ-l.NPM 

IF(YP(I J).LE.10.A*0NHIN) YPCIJ) - 10. ♦♦DNMIN 
YPIIJ) - ALOGIO ( YP ( I J ) ) 

600 CONTINUE 

IF( IPSLCT(6).NE.O) 

1CALL AXE S (0. *0. *0.*8.«X(NP*l)*X(NP+2) *1.*0. * 

1 1HX * 0. 2 *-l ) 

IF ( I PSLCTI 6 ) . NE . 0 ) 

1C ALL AXES ( 0. *0.*90. *6. «YP ( NP ) * YP ( HP* l ) ,Z.*0. . 

I 6HL0G { J I *0.2,6) 

IF(IPSLCT(15).NE.0) 

1CALL AXESC-2.*0. *90.*6.,YP(NP) *YP INP+1) ,2.* 0 . 9 

1 *N ( UCLE AT ION )R { ATE » S2LS30G )J* D ( ROP LETS /MS U3SN-S ) 1 * 

2 0.2,51) 

CALL LI NPLT ( X( 2) *YP* NPM, 1*0* 0*0*0) 

CALL NFRAME 

CALL CALPLT (A.5,1.5,-3) 

610 IFdPSLCT(ll).EO.O) GO TO 630 
C 

C ENTROPY PRODUCTION VS X 

C 

DO 620 IJ-1,NP 
YPCIJ) - SIIJ) - SO 
620 CONTINUE 

CALL ASCALE(YP*6.*NP«1 *10. ) 

IF { SMIN.NE.0..0R.5MAX.NE.0.) YP<NP*1) - SMIN 
IFISMIN.NE.O..OR.SMAX.NE.O.) YP(NP*2) - (SMAX - SMIN)/6. 
I F ( IPSLCT (ll).NE.O) 

1C ALL AXES(O.*0.*O.*8.*X(NP*1) , X (NP + 2 ) *1 • * 0 . * 

2 'DUSTANCE ) A (LONG )N(OZZLE X, M)', 0.2, -321 
IF C IPSLCT (11) . NE • 0 ) 

1C ALL AXES (0. ,0. ,90. *6. ,YP(NP*1) *YP(NP*2) *1. *0. * 

1 'EiNTROPY >P( ROOUCTION, )SA043(S, )J/|KG)-K»* 

2 0 .2 * A 2 ) 
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503 

504 

505 

506 

507 
500 

509 

510 

511 

512 

513 

514 

515 

516 

517 
510 
519 


CALL LINPLT(X,YP,NP,1»0,0»0,Q} 

CALL NF RARE 

CALL CALPLT(4.5,1.5,-3> 

CALL N0TATE{0**5.5*0.16t , TSD(S) SN($(PS01»NS)-m0GSN' »0.,27) 

CALL NOTATE (0. *4.5,0. 16 »« J'tO. ,1) 

CALL NO TATEIO. *3.5,0.16* • <Ml t i0.i3) 

CALL NOTATE (0. ,2.5*0.16, 1 ( SATURATION ) 1 ,0.*12) 

CALL NOTATEtO. ,1.5,0.16, * ( R ) tDMN 1 , 0 . * 8 ) 

CALL LINPLTC PATX.PATY ,2, 1*0, 0,0*0) 

CALL N0TATE{1.7,O.5f0.16,' (NUCLE AT ION AND CONDENSATION >■, 0. * 29) 

PATY(1I - PATYUJ - 0.5 

PATYiZl - P A TY ( 2 ) - 0.5 

CALL DASHPLT(PATX,PATY,2,1,PAT,2) 

CALL NOTATEtl. 7,0. ,0.16, MI5ENTR0P1C WITHOUT CONDENSATION I ■ ,0 ., 33 ) 
630 CONTINUE 
RETURN 
END 
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SUBROUTINE PRTQUT ( 1 P » IL ) 

C 

C SUBROUTINE PRTOUT GENERATES PRINT OUTPUT OF THE 

C PROGRAM RESULTS 

C 

DIMENSION HLI3 (400) ,$LI3 (400) ,YI J3 (400) 

COM MON/CM 0UT3/NENT 
COKHON/GAMGUT/GAMBAR 

COMMON/ INP T/TO , PO.XLASTt PH II, RHC*RNO 

COMMON I3,J3,A3(400) ,X3 (400) ,DELN3(400) , CAP JX3 ( 400 ) , R3 (400) , 

1 MU3(400) ,TL3(400) ,T3(400) , P3 ( 400 ) , RST AR3 ( 4 00 ) , 

2 U3(400 ) ,S3 (400),SL3(400) ,HL3(400) ,RH03(400) , 

3 CAPJY3(400) ,F3( 400) ,G3( 400) ,YS 3(400 ) , 

4 DELY3( 400) , DELS 3 (400) ,RPERM3(400) ,TS3(400) 

5 ,OUTMCH( 400) ,CAPP 3(400) 

REAL MU3 

COMMON /OUT 1/ SO ,H0 , MO OT , RHO S TR , T S TAR 
REAL MOOT 

COMMON/ CMOUTl/NENTRYf DT*RHOO f Tl#HltU It US TAR ,URHO » RHOF 1 
COMMON/S IGC0M/SIGMA(10),SIG,PSUP0 

COMMON/ 0UT2/U2»RHO2*CAPJ2,T2,P2,RSTAR2,X2,A2,SIG2(il), OUMC 
COMHON/AOJUST/GAMMAI »QC, ALPHA, BETA, DSTAR, AS TAR # DT2, 

1 JMIN,DELX,EPS(8),I0UT,DT1,CVRAT,ARAT, JSTOP , JOB,CP ,CPl 
REAL JM IN 

C0MM0N/SWITCH/SWEND,SWP0,5*GQ 

C0MMDN/CM0UT2/NE,T2«,SIGMAX,U2M,RH0,P,PINF,PBAR,TR,H2N,RSTAR,CAPL, 

1 PH I , J 
REAL J 

INTEGER SWENDtSWPOiSwGO 
DATA N2/0/ 

DATA IFLG/1/ 

C 

1 FORMAT (1H1,20X,14HPART 1 SUMMARY/// 

1 1H0 1 17X , 9HS 0 « , 1PE 16 *8 / 

2 1H0,17X,9HH0 - .1PE16.8/ 

3 1H0, 17X,9HMD0T - ,1PE16.8/ 

4 IHO , 1 7X , 9HRH0* * ,1PE16.8/ 

5 IHO * 1 7X , 91-T* « ,1PE16.8/ 

6 IHO , 17X , 9HU* - ,1PE16.8/ 

7 1H0,17X,9HS1GMA - ,1PE16.8/ 

8 9(1H0,26X,1PE16.8/) ) 

2 FORMAT (1H1 ,20X ,22HPART 1 EXTENDED OUTPUT/// 

1 1 HO, 4 X , 4X , 5H A* ,9X,3HH0 ,9X,5HRH00 ,9X,3HS0 ,10X, 

2 3HRH0, 1 IX , 1H T , 12X, 1HH , 12 X , 1HU , 10X , 5 HU* RHO , 8X, 6HGAMMA1/ ) 

3 FORMAT (3X,10(1PE13*5 ) ) 

4 F0RHAT(1H1,20X,14HPART 2 SUMMARY/// 

1 1H0,17X,9HU - , 1PE 16*8/ 

2 IHO, 17X,9HRH0 « ,1PE16.8/ 

3 1H0,17X,9HJ - , IP E 16 * 8/ 

4 IHO , 1 7X , 9hT - » 1PE16 • 8/ 

5 IHO , 17X * 9HP - , 1 PE 16 • 8/ 

6 IHO , 17X ,9HR4 - ,1PE16.8/ 

7 IHO, 17X.9H51CMA - ,1PE16.B/ 

8 9(1H0,26X,1PE16.8/) ) 

5 FORMAT (1H1,20X,22HPART 2 EXTENDED OUTPUT//) 

6 FORMAT ( _ 

1 1H0,4X,4X,4HT/T0,11X,1HH,12X,1HU,11X, 3HRH0 ,9X, 4HP / PO ,9X,5H PINF, 

2 8X,4HPBAR,11X, 1HX , 10X , 6HG AMMA1 ) 

7 FORMAT ( 

1 1H ,4X,5X*3H TR,9X,6HSIGMAX,8X,5HRSTAR,10X,1HL,11X,3HPHI , 

2 11X,1HJ) 

8 FORMAT ( 1H1,2CX*14HPART 3 SUMMARY/// 

1 1H0,4X,6X,1HX,1QX,4HP/P0,9X,4HT/T0, 

2 1QX,2HJ-,11X,2HJ+»9X , 4HDELN , 1 IX • 1HU , 12X ,1HY, 10X,4HDELS/ ) 

9 FORMAT (1H1,20X,22HPART 3 EXTENDED OUTPUT//) 

10 F0RMAT(I4,9(1PE13.5>) 

11 FORMAT ( 1H1 ) 

12 F0RMAT(//5X,4HJ - , I 5, 5X , 4HX * ,1PE13*5, 

2 5X,4HT * ,1PE16*8//7X,2X, 1HI, 

1 8X , 1HR , 1 IX , 2HMU , 1 IX , 3HHLI , 10X,3H$LI , 10X.2HTL, 10X,3HYI J/) 

13 F0RMAT(5X,I5,6(1PE13.4J) 

14 FORMAT ( 5H R - , 1 X ,9 ( IP E 1 3 . 5 ) /30 ( 6X , 9 ( IP El 3 • 5 ) / ) ) 

15 FORMAT! 

1 1H ,2X*1HJ»1X,3X»6HT( J* 1 ) , 7X, 6HP (.)♦ 1) ,7X,6HUU + 1) , 7X, 6HY ( J* 1 ) , 

2 7X,6HS(J*1) ,7X,7HSL(J*1) ,6X, 7HHL ( J* 1 ) ,6X,6HX( J*l) »7X, 6HA ( J *1) ) 

16 F0RMAT(5X,2X,8HRH0(J*l) , 4X, 10HR5T ARC J*1),4X,9H DELYU) , 

1 4X,9H DELN(J) ,5X,7HTS ( J*l) »4X, 8HCPH J*l> , 7X,6HGAMHA1) 

17 F0RHAT(4X,9(1PE13*5) ) 

C 

IFCNENTRY.EG.l) NLINE - 0 
GO TO (20,40,70) ,IP 
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83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 


C PART 1 OUTPUT 

C 

20 CONTINUE 

IF ( IL .ME • 0) GO TO 30 
C 

C STANDARD OUTPUT 

C 

WRITE (6,1) SO, HO, HDOT, RHOSTR.TS TAR, U STAR, SIGMA 
GO TO 140 
30 CONTINUE 

IF1IL.NE.2) GO TO 140 
C 

C DEBUG OUTPUT 

C 

IF(NENTRY.EQ.l.OR.NLINE.EQ.O) WRITE (6,2) 

IF(NLINE.EQ.C) NllNE - 6 

NLINE - NLINE ♦ 1 

IF ( NLI NE • GT • 45 ) NLINE - 0 

WRITE (6,3) ASTAR,H0,RH00,S0,RH0F1,T1,H1,U1,URH0,CAMBAR 
GO TO 140 

c 

C PART 2 OUTPUT 

C 

40 CONTINUE 

IF(IL.NE.O) GO TO 50 
C 

C STANDARD OUTPUT 

C 

WRITE (6, 4) U2* RH02 »CAP J2 , T2 , P2 , RSTAR2 , S IG2 
GO TO 140 
50 CONTINUE 

IF(IL.NE.l) GO TO 60 
C 

C EXTENDED OUTPUT - 1 

C 

IF(NENTRY.E0.1 .OR. NLINE .EQ.O) WRITE (6, 5) 

I F ( NENTR Y *L T. N2 ) I FLG - 2 
IF(NLINE.EQ.O) WR I TE ( 6 i 6 ) 

IF (NLINE .EQ.C. AND. IF LG .EG .2) WRITE (6, 7) 

IF(NLINE.EQ.O) NLINE - NLINE ♦ 6 
IFCNLINE.EQ.6.AND.IFLG.EQ.2) NLINE - NLINE ♦ 1 
NLINE - NLINE ♦ 1 


125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 


IF (NLINE.GT • 4 5 • AND • IFLG.EQ.l) NLINE » 0 
XTT - T2W/T0 
PP - P/P 0 

WRITE (6, 3) XTT,H2W,U2W,RHO,PP,PINF,P0AR,X2, GAM BAR 
N2 - NENTRY 
GO TO 140 
60 CONTINUE 
C 

C EXTENDED OUTPUT - 2 

C 

NLINE • NLINE ♦ 1 
IF ( NLI NE • GT *45) NLINE - 0 
WRITE (6,3) TR,SIGMAX,RSTAR,CAPL,PHI, J 
GO TO 140 
C 

C PART3 OUTPUT 

C 

70 CONTINUE 

IFCIL.NE.O) GO TO 120 
C 

C STANDARD OUTPUT 

C 

WRITE(6»8) 

NLINE - 6 
J3M1 » J3 - 1 
DO 80 I-1.J3M1 
IF(NLINE.EQ.C) WRITE(6,8) 

I F (NLINE . EO • 0 ) NLINE - NLINE ♦ 6 

NLINE - NLINE ♦ 1 

IF (NLI NE.GT.45) NLINE - 0 

PP - P3(I)/PG 

TT - T3 ( I ) /TO 

WRITE (6,3) X 3 ( I ) , PP , TT, CAP JX3 (I ) , CAP JY3 ( I.) ,DELN3 ( I ) , U3{ I ) , YS3 ( I) 
1 .DELS 3(1) 

00 CONTINUE 

I F ( IL.NE.IOO) GO TO 120 
REWIND 7 
WRITE (6, 11) 

90 CONTINUE 

READ! 7, END- 100) J J J ♦ ( R3 ( I ) ,MU3 ( I ) , HL 1 3 ( I ) , SL I 3 ( I ) , TL 3 ( I ) , YI J 3 ( I ) 

1 , I -1 , J J J ) 

100 IF (EOF (7) .NE • 0) GO TO 140 
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167 

166 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
162 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 


*RITE(6,12> JJJ*X3( JJJ) .T3UJJ) 

do no i-i,jjj 

WRITE* 6 , 13) I«R3III»MJ3in *HL 1 3 ( I > * SL I 3 ( I > * TL 3 ( I ) * YI J3(I J 
110 CONTINUE 
CO TO 90 
120 CONTINUE 
C 

c DEBUG OUTPUT 

C 

IFUL.NE.l) GO TO 130 
IFU3.EQ.il NUNE * 0 
IF ( NLI NE.EQ . 0) WRITE(6,9) 

IF(NLINE.EQ.O) WRITE(fe*15) 

IF(NLINE.EQ.O) WRITE<6,16) 

IFtNLINE.EQ.O) NLINE - NLINE ♦ 9 ♦ (J3^1)/9 

WRITE (6 *10) J3,T3< J3 + 1 I t P 3 ( J3*l ) *U3( J3*l) ,YS3( J3*i) »$3(J3*1) » 
1 SL3 ( J 3 + 1) *HL3 ( J3 + 1 ) * X3 < J3 ♦ 1 J t A3 ( J3* 1 1 
WRITE (6, 17) RH03 ( J3+ 1) *RSTAR3( J3+1I*DELY3(J3)*DELN3( J3>, 

1 T$3(J3*1)*CP1*GAMBAR 
WRITE {6 *14) ( R 3 ( I JKL)* I JKL-1* J3) 

NLINE - NLINE ♦ 3 ♦ (J3*l)/9 
IF (NLINE. GT. 45) NLINE - 0 
GO TO 140 
130 CONTINUE 
140 CONTINUE 
RETURN 
END 


1 

2 

3 

4 

5 

6 
7 
6 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 


SUBROUTINE S AT T EM ( T5 * PX ) 

C 

C SUBROUTINE SATTEM CALCULATES THE SATURATION TEMPERATURE 

C 

COMMON/S IGC0M/SIGMAU0),$IG,P5UP0 

COMMON/ AD JUS T/GAMMA1, QC, ALPHA, BETA, DSTAR, AS TAR, DT2i 
1 JMIN*DELX,EPS(8), I OUT *DTl*CVRAT*ARATtJ$TOP,JDB,CPiCPl 
REAL JM IN 
C 
C 

C SOLVE EQUATION (57) FOR T 

C 

TOLD - 300. 

ALOP - AL0G(PX*SIGMA(1)/SIG) 

10 CONTINUE 

F ■ 55.897 - 6641.7/TOLD - 4. 48 64*ALQG ( TOLD ) - ALOP 
FP - 664 1.7/T0LD4*2 - 4.4864/T0L0 
TNEW - TOLD - F/FP 

IF(ABS(TNEW-T0LD)/T0LD.LT.EPS(1)) GO TO 20 
TOLD - TNEW 
GO TO 10 
20 CONTINUE 
TS - TNEW 
RETURN 
END 
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ORIGINAL PAGE IS 
OF POOR QUALITY 

1 SUBROUTINE XVSA < ALST .XLST ) 

2 C 

3 C SUBROUTINE XVSA CALCULATES THE X POSITION IN THE 

A C NOZZLE AS A FUNCTION OF AREA 

5 C ACOEF CONTAINS THE SPLINE COEFFICIENTS DETERMINED 

b C BY THE LIBRARY ROUTINE CSOS 

7 C 

8 COMMON/ARECOM/NA,XA,ACOEF 

9 DIMENSION ACOEFtlOfAJ 


10 


DIMENSION A I 10 ) 

11 


COMMON/ CONST /TC,K,CAPR,M1,W( 10 ) ,RH0L #RBAR f P IfALPHAC 

12 


REAL K,M1 

13 


COMMON/ AD JU ST/ 6 AMMAltQCf ALPHA t BETA fDSTAR 9 AS TAR fDT 2 « 

1 A 


1 JMIN, DELXf EPS( 8) 1 1 OUT , DTI * C VR AT t ARA T » JSTOP *JDBtCPfi 

15 


REAL JM IN 

lb 


DIMENSION XA(10) 

17 

C 


18 


DO 10 1*1, NA 

19 


A(I) - AC OEF ( 1,1) 

20 


MI) - PI*AIX)**2 

21 

10 

CONTINUE 

22 


AREA - ALST 

23 


RBARX • SORT (AREA/PI) 

2 A 


NPD - NA-1 

25 


DO 20 I J* 1 , NPD 

26 


IFUREA.GE.AU J). AND. AREA. LE.AUJ + l)) GO TO 30 

27 

20 

CONTINUE 

28 


IJ - NPD 

29 

30 

CONTINUE 

30 


HOLD • 0.5*1 XAUJ+1) - XA(IJ)) 

31 

AO 

CONTINUE 

32 


FXX - t(ACOEF(IJ,A)*HOLD ♦ ACOE F { I J , 3 ) ) AHOLD 

33 


1 ♦ ACOEFt I J,2 ) )*HOLD ♦ ACOEFUJtl) 

3 A 


FXX - FXX - RBARX 

35 


FXXP - C3.*AC0EF(I J,4)*H0LD ♦ 2 , * ACOE F ( IJ , 3 ) ) *HOLD 

36 


1 ♦ AC0EF(IJ,2) 

37 


HNEW - HOLD - FXX/FXXP 

38 


IFUBS (HNEW-H0LD)/H0LD.LT.EPS(8) ) GO TO 50 

39 


HOLD - HNEW 

AO 


GO TO 40 

A 1 

50 

CONTINUE 

A2 


XLST - XA( I J) ♦ HNEW 

A3 


IFUREA.LT. Ail) ) XLST « ( AREA/ A( 1) ) *XA (1 ) 

AA 


RETURN 

A5 


END 


61 



References 

1. Hill, Philip G.: Condensation of Water Vapor During 
Supersonic Expansion in Nozzles. J. Fluid Mech ., vol. 25, 
pt. 3, July 1966, pp. 593-620. 

2. Gyarmathy, G.: Condensation in Flowing Steam. Two- 
Phase Steam Flow in Turbines and Separators, M. J. 
Moore and C. H. Sieverding, eds., Hemisphere Publ. 
Corp., c.1976, pp. 127-189. 

3. Young, J. B.: Spontaneous Condensation of Steam in 
Supersonic Nozzles. Dep. of Engineering, Univ. of Cam- 
bridge, 1980. 

Part I — Nucleation and Droplet Growth Theory. 
CUED/A-Turbo/TR 97, Pt. I. 

Part II — Numerical Methods and Comparison With 
Experimental Results. CUED/A-Turbo/TR 97, Pt. II. 

4. Young, J. B.: The Spontaneous Condensation of Steam 
in Supersonic Nozzles. PhysicoChem. Hydrodyn ., vol. 3, 
no. 1, 1982, pp. 57-82. 


5. Keenan, Joseph H.; Keyes, Frederick G.; Hill, Philip 
G.; and Moore, Joan G.: Steam Tables — Thermodynamic 
Properties of Water Including Vapor, Liquid, and Solid 
Phases. John Wiley &; Sons, Inc., c.1978. 

6. Kantrowitz, Arthur: Nucleation in Very Rapid Vapor 
Expansions. J . Chem. Phys., vol. 19, no. 9, Sept. 1951, 
pp. 1097-1100. 

7. Touloukian, Y. S.; Saxena, S. C.; and Hestermans, P.: 
Viscosity. Thermophysical Properties of Matter, Vol- 
ume II, IFI/Plenum, c.1975. 

8. Touloukian, Y. S.; Liley, P. E.; and Saxena, S. C.: 
Thermal Conductivity — Nonmetallic Liquids and Gases. 
Thermophysical Properties of Matter, Volume S, 
IFI/Plenum, 1970. 

9. Anderson, E. C.; and Lewis, C. H.: Laminar or Turbu- 
lent Boundary- Layer Flows of Perfect Gases or Reacting 
Gas Mixtures in Chemical Equilibrium. NASA CR-1893, 
1971. 

10. JANAF Thermochemical Tables. Thermal Lab., Dow 
Chemical Co., June 30, 1965. 


62 



Table 1. Conditions for Numerical Results 
[For all cases: q c = 1, a = 8, a c — 1, /? — 2, and J m i n — 10 15 ] 


Case 

To, K 

Po, bars 

<t> 

yH 2 0 

1 

1900 

50 

0.798 

0.154 

2 

1900 

250 

.797 

.154 

3 

1600 

50 

.620 

.122 

4 

1600 

250 

.620 

.122 

5 

1900 

250 

.427 

.156 


Table 2. Nozzle Coordinates for Langley 8'HTT and Computed 
Boundary-Layer-Displacement Thickness 


x, m 

r Wy m 


6* y m ) for- 

— 


Case 1 

Cases 2 and 5 

Case 3 

Case 4 

0 

0.0710 

0 

0 

0 

0 

.0860 

.0714 

0 

0 

0 

0 

2.5400 

.3538 

.0072 

.0059 

.0073 

.0060 

5.4703 

.6911 

.0271 

.0213 

.0275 

.0217 

6.4230 

.7938 

.0358 

.0281 

.0363 

.0286 

9.8560 

1.0434 

.0670 

.0523 

.0675 

.0528 

12.8949 

1.1660 

.0911 

.0713 

.0916 

.0718 

15.7116 

1.2192 

.1089 

.0851 

.1093 

.0855 








Normalized Static Pressure, p/p 0 Nucleation Rate, fog J, Droplets/m 



O 2 4 6 8 10 12 14 16 

Distance Along Nozzle x, m 


(a) Temperature difference, nucleation rate, 
and mass fraction liquid water. 


12 x 10' 4 



0 2 4 6 8 10 12 14 16 


Distance Along Nozzle x, m 

(c) Static pressure distribution. 

Figure 2. Numerical results for case 1 with T 0 


Entropy Production, As, J/kg-K 


30 



(b) Entropy production. 



0 2 4 6 8 10 12 14 16 

Distance Along Nozzle x, m 


(d) Droplet growth. 

1900 K, p 0 = 50 bars, and <fi = 0.798. 
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Temperature Difference, T^p^-T 




Distance Along Nozzle x, m 


(a) Temperature difference, nucleation rate, 
and mass fraction liquid water. 


(b) Entropy production. 




(c) Static pressure distribution. 


(d) Droplet growth. 


Figure 3. Numerical results for case 2 with T 0 = 1900 K, p Q = 250 bars, and (j> = 0.797. 


Normalized Static Pressure, p/p Q Nucleation Rate, fog J, Droplets/m 3 -: 



(a) Temperature difference, nucleation rate, 
and mass fraction liquid water. 


O' 25 
\ 


Q_ 10 



J I h _ _1 I J 

2 4 6 8 10 12 

Distance Along Nozzle x, m 


L 


14 16 


(b) Entropy production. 




(c) Static pressure distribution. (d) Droplet growth. 

Figure 4. Numerical results for case 3 with T 0 = 1600 K, p Q — 50 bars, and <f) = 0.620. 
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Normalized Static Pressure, p/p G Nucleation Rate, tog J, Droplets/m 



O 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 

Distance Along Nozzle x, m Distance Along Nozzle x, m 


(a) Temperature difference, nucleation rate, (b) Entropy production, 

and mass fraction liquid water. 



0 " " 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 

Distance Along Nozzle x, m Distance Along Nozzle x, m 


(c) Static pressure distribution. (d) Droplet growth. 

Figure 5. Numerical results for case 4 with T 0 = 1600 K, p 0 = 250 bars, and <j> = 0.620. 
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Normalized Static Pressure, p/p 0 Nucleation Rate, tog J, Droplets/m 



(a) Temperature difference, nucleation rate, 
and mass fraction liquid water. 



(b) Entropy production. 




10b 


8b 
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\ 

\ 

\ 


nucleation and condensation 

isentropic without condensation 


saturation 



o I 1 l i i J I i J 

O 2 4 6 8 10 12 14 16 

Distance Along Nozzle x, m 



(c) Static pressure distribution. 


(d) Droplet growth. 


Figure 6. Numerical results for case 5 with T 0 = 1900 K, p 0 = 250 bars, and 4> — 0.427. 
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